NASA-CR-196292 




/ S3 o 4 > 


FINAL REPORT 


NASA CONTRACT NAG 8-149 
PROJECT JOVE 


AUGUST, 1994 


(NASA-CR-196292) PROJECT JOVE N94-37635 

Final Report (West Virginia Uni v. ) 

116 p 

Unci as 


G 3/34 0018306 


Dr. M.J. Lyell, p.i. 
Mechanical & Aerospace Engr. Dept. 
Box 6106 

West Virginia University 
Morgantown, WV 26505-6106 
(304) 293-3111 



FINAL REPORT 


NASA CONTRACT NAG 8-149 
PROJECT JOVE 


AUGUST, 1994 


Dr. M.J. Lyell, P.I. 
Mechanical & Aerospace Engr. Dept. 
Box 6106 

West Virginia University 
Morgantown, WV 26505-6106 
(304) 293-3111 



CONTENTS -• 


Title page 1 
Contents 2 
Project overview 3 
Project Accomplishments 7 


Section A: Investigation of Potential Research Topics of 
Scientific and Commercial Interest Relevant to 
Microgravity Fluid Mechanics 


Section B: Research Results: 

Interface Response of the Finite Length Liquid 
Column to Longitudinal Forcing 
and 

Interface Stability of the Infinite Liquid 
Column to Forcing Normal to the Column Axis 


Section C: Finite Length Liquid Column Stability 
in the Presence of a Periodic 
Acceleration Field Oriented Normal 
the Longitudinal Axis 


Section D: Copy of the Publications made under Project 
JOVE, NASA Contract NAG8-149 

Section E: Copy of Syllabus — Educational Effort 


Supplement I: Interface Behavior of a Multi-Layer Configuration 

Subject to Acceleration in a Microgravity 
Environment 

(THIS IS IN A SECOND SPIRAL BOUND VOLUME) 


Supplement II: Nonlinear Effects on the Natural Modes 

of Oscillation of a Finite Length 
Inviscid Fluid Column 
(THIS IS IN A THIRD SPIRAL BOUND VOLUME) 


2 



PROJECT OVERVIEW: 


This final report for Project JOVE (NASA Contract NAG8-149) 
presents the results of work performed. Efforts under this 
contract involved a research component, an education component, 
and an outreach component. 

Early on, the research effort under JOVE was concerned with 
identifying specific research activities of interest to both the 
principal investigator and NASA technical interests. In the work 
of this principal investigator, the research efforts were 
concerned with the behavior of fluid dynamics in a microgravity 
environment; in particular, the free surface/ interface behavior of 
fluid configurations. Earlier work centered on the behavior of 
"slab” configurations. More recent efforts have looked at the 
finite length fluid column behavior under forcing. It is 
emphasized that the multiplicity of research tasks was considered 
to be open-ended (ie., there's always more research to be done). 

Indeed, part of the agenda of the JOVE project has been to 
involve more university researchers in research of interest to 
NASA. To this end, it was understood that the principal 
investigator would seek additional sources of funding. The 
proposals which were written to do this are listed in the 
following chapter on Project Accomplishments. 

A new focus for Year 4 of the JOVE project was the effort to 
identify practical applications of free surface fluid mechanics 
which would have applications in microgravity fluid mechanics or 
be relevant to another aspect of the NASA mission and be also of 
some commercial interest. The thrust was to seek not only novel 
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application area(s) in which research would contribute to the 
technology, but aiso to identify problems which are accessible via 
a fully computational fluid dynamics methodology. This is a 
departure from the previous JOVE efforts, which involved primarily 
theoretical formulations. The results of this investigation are 
listed in Section A of this volume of the final report. 

Research results on the interface shape and stability of the 
slab configuration are given in the additional volume labeled 

"SUPPLEMENT I INTERFACE BEHAVIOR OF A MULTI-LAYER FLUID 

CONFIGURATION SUBJECT TO ACCELERATION IN A MICROGRAVITY 
ENVIRONMENT". (This volume is essentially the thesis work done by 
a graduate student who was supported by the contract, and who 
received a masters degree in mechanical engineering.) 

With regard to the research investigations which focused on 
the interface behavior of the fluid cylinder in microgravity, the 
results may be perused as follows. The problem of the response of 
a finite length fluid column to forcing along the longitudinal 
axis of the column was accomplished via Laplace transform 
techniques. This work is presented in Section B of this volume of 
the final report. The results of the investigation on the 
interface stability of an infinite length fluid column subject to 
periodic accelerations oriented perpendicular to the longitudinal 
axis of the column are found in Section B also. 

Nonlinear behavior was taken into account in one of the 
investigations of finite length fluid column interface behavior in 
microgravity. In this case, the focus was on the nonlinear 
corrections to the normal modes of free oscillation of the finite 
length fluid column. It was not necessary that the fluid column be 
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restricted to the slender column limit case. The analysis was 
inviscid and incompressible, and the perturbations were assumed 
irrotational. Also, the disturbances were restricted to those 
which are axisymmetric . Additional details and results on this 

problem are found in "SUPPLEMENT II NONLINEAR EFFECTS ON THE 

NATURAL MODES OF OSCILLATION OF A FINITE LENGTH INVISCID FLUID 
COLUMN". This is in another spiral volume of the final report. 
(This volume is essentially the first draft of the dissertation of 
a Ph.d. student who has been supported by the JOVE contract.) 

With regard to the problem of the interface stability of the. 
finite length fluid column subject to both periodic disturbances 
and a static acceleration oriented normal to the longitudinal axis 
of the column, the full theoretical formulation of the problem, 
including the static deformation, has been made. The numerical 
investigation is underway, and when finished, results will be 
submitted for publication. A discussion may be found in Section c 
of this volume of the final report. 

The educational efforts of this program have involved the 
instruction and supervision of graduate students ( 2 masters level 
and 1 ph.d level), and the teaching of a graduate course 
specialized to include free surface/ interface behavior in 
microgravity. These will be listed in the JOVE accomplishments. A 
syllabus for the graduate class is attached in Section E of this 
volume of the final report. 

The outreach efforts of this project consisted primarily of 
lectures/talks given to the general public. However, one of the 
outreach events involved a "workshop style" presentation to middle 
school- junior high age girls. The workshop introduced the girls to 


5 



the world of an aerospace engineer, and involved some hands-on 
demonstrations of flight principles. Spacecraft and microgravity 
were part of this. (The workshop was co-sponsored by the WV 
chapter of the Association for Women in Science (AWIS) , the WV 
Dept, of Education — Office of Sex Equity, the AAUW, and the 
Southern West Virginia Community College) . 
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JOVE ACCOMPLISHMENTS 


RESEARCH: 


REFEREED JOURNAL ARTICLES PUBLISHED: 

(1.) M.J. Lyell, "Fluid column stability in the presence of 

periodic accelerations", AIAA Journal, Vol. 31, p. 1519-21, 1993. 

(2.) M.J. Lyell, "Axial forcing of a finite length liquid 

column", Phys. Fluids , Vol. 3, p. 1828-31, 1991. 

(3.) M.J. Lyell and M. Roh, "The effect of time-dependent 

accelerations on interface stability in a multi-layered 
configuration", AIAA Journal, Vol. 29, p. 1894-1990, 1991. 

Manuscript in Preparation; to be submitted to Physics of Fluids 
by Sept. 30, 1994, entitled, "Nonlinear corrections to interface 

shape and oscillation frequencies of a finite length inviscid 
liquid column in microgravity", by M.J. Lyell and L. Zhang. 

OTHER PUBLICATIONS: 

(1.) M.J. Lyell and L. Zhang, "On the nonlinear dynamics of liquid 
bridges", Proceedings, 12th U.S. National Congress of Applied 
Mechanics, 1994. (Conference held in Seattle, WA, in June, 1994.) 
(2.) M.J. Lyell, "Interface stability of a fluid column subject to 
periodic acceleration oriented normal to the longitudinal axis of 
the column", presented at the World Space Congress, held 
Washington, D.C., Aug. 1992, PAPER NUMBER IAF-92-0914. 

(3.) M.J. Lyell, "Interface stability of a liquid column in the 
presence of periodic accelerations oriented normal to the 
longitudinal axis", Bull. Am. Phys. Soc., Vol. 37, 1992. 
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(4.) M.J. Lyell and M. Roh, "Instabilities in a multi-slab fluid 
configuration due to time-dependent acceleration normal to the 
fluid-fluid interface in a microgravity environment", 29th AIAA 
Aerospace Science Conference (held Jan., 1991, Reno, NV) , 

AIAA PAPER 91-1019. 

(5.) M. Roh and M.J. Lyell, "Investigation of Interface in a 
Multiple Layer Slab Configuration — Utilizing WVNET Computational 
Resources", WVNET Conference 1990 Proceedings, p. 27-39, 1990. 
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TEACHING: 


GRADUATE STUDENTS : 

(1.) Ms. L. Zhang 

Ph.d. Candidate, MAE 

Research Area: fluid mechanics/non-linear oscillations. 
Dissertation Topic: Nonlinear corrections to the natural 

oscillations of a finite length inviscid liquid column in 
microgravity. 

Anticipated graduation date: Dec. 1994. 

(2.) Ms. K. Perkins 
MSAE Candidate 

Research Area: Fluid mechanics , ferrofluids, free surfaces. 
Thesis topic: Wave dynamics in ferrofluids. 

Anticipated degree date: Dec. 1995. 

(3.) Mr. Michael Roh 

MSME Degree awarded May, 1991. 

Research area: Fluid mechanics. 

Thesis topic: Stability of fluid layer configurations subject 
to time-varying acceleration with application to microgravity 
fluid mechanics. 

COURSES: 

(1.) MAE 399 —SPECIAL TOPICS GRADUATE COURSE 

entitled "Fluid dynamics of free surfaces/ interfacial 
fluid mechanics" . 

Course formally taught to 1 student in Fall semester, 1993. 
Syllabus for this course in Section E of this volume of the 
final report. 
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OUTREACH : 


( 1 . ) Event : 


(2.) Event: 


( 3 . ) Event : 


Workshop and Question and Answer Panel. 

March, 1994 , at WV Southern Community College, 
Williams, WV. 

65-80 participants (middle school- junior high girls) . 
workshop on "Up, up, and Away, The World of an 
Aerospace Engineer" . 

Lecture/ film / handouts 

Oct., 1991, at Cheat Lake Junior HS, Morgantown (Cheat 
Lake area) , WV. 

80 participants. 

handout on space suit designs. 

Invited speaker at Sigma Xi meeting. 

Dec., 1990, at Marshall Univ. , Huntington, WV. 

25-30 participants. 

"general lecture" to audience with diverse scientific 
backgrounds . 
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SECTION A 



Investigation of Potential Research Topics of 
Scientific and Commercial Interest Relevant to 
Microgravity Fluid Mechanics 



Abstract 


The goal of this project is to investigate new areas of 
research pertaining to free surface-interface fluids mechanics 
and/or microgravity which have potential commercial applications. 
This paper presents an introduction to ferrohydrodynamics (FHD) , 
and discusses some applications. Also, computational methods for 
solving free surface flow problems are presented in detail. Both 
have diverse applications in industry and in microgravity fluids' 
applications. Three different modeling schemes for FHD flows are 
addressed and the governing equations, including Maxwell's 
equations are introduced. In the area of computational modeling of 
free surface flows, both Eulerian and Lagrangian schemes are 
discussed. The state of the art in computational methods applied 
to free surface flows is elucidated. In particular, adaptive grids 
and rezoning methods are discussed. 
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Chapter I 
Overview 


The goal of this project was to investigate new areas 
pertaining to microgravity which have potential commercial 
applications. The two topics studied for this project were state 
of the art computational methods for solving free surface flow 
problems and ferrohydrodynamic free surface flows. 

Free surface and interface flows do indeed have many 
commercial applications in diverse areas. These include crystal 
growth, melting and solidification, capillary flows, wave 
propagation, and metal and glass forming processes. 

Ferrohydrodynamics deals with fluid dynamics that occurs in a 
magnetic fluid as a result of an applied magnetic field. A 
ferrofluid is a colloidal suspension of solid magnetic particles in 
a typically Newtonian parent liquid. These fluids also have many 
applications, including rotary shaft seals, cooling processes, and 
sink-float separation processes. 

A number of interesting phenomena are exhibited by magnetic 
fluids in response to applied magnetic fields. These include a 
normal field instability, in which a pattern of spikes appears on 
the fluid surface. Also, in thin layers of ferrofluid, there 
exists the spontaneous formation of labyrinthine patterns. In 
rotary fields, a body couple is generated which is manifested as 
anti-symmetric stress. In addition, ferrofluids with temperature- 
dependent magnetic moment allow for enhanced convective cooling. 
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Ferrohydrodynamics is a fairly new area of study, since the 
onset of research in this field began in the 1960s. There is still 
much to be discovered in this area, and linking ferrofluids with 
free surface flow problems can lead to many interesting new 
developments . 

Chapter II of this paper gives a brief introduction into 
ferrohydrodynamics. This section also discusses some of the many 
applications involving ferrofluids. 

There are three primary modeling schemes which are appropriate 
in analyzing ferrofluid motion — (1) restriction to inviscid flow, 
(2) viscous flow with a symmetric stress tensor, and (3) viscous 
flow with a non- symmetric stress tensor. Situations are given in 
which each scheme would be applicable. In addition, the main 
governing equations for ferrofluid motion in each modeling type and 
the equations for the magnetic field conditions (including 
Maxwell's equations) are included in the discussion. 

Chapter III discusses state of the art computational methods 
for solving free surface flow problems. Found in this chapter is 
a discussion on Eulerian and Lagrangian computational methods for 
solving free surface problems, with several exairqples given in 
detail. Also discussed is the efficacy of finite differences 
versus finite element formulations. 

Chapter IV gives conclusions regarding this project. 
Suggestions are provided for future research ideas arising from the 
investigation of free surface flows and ferrohydrodynamics. 



Chapter II 

Introduction to Ferrohydrodynainics 


Ferrohydrodynamics (FHD) deals with the mechanics of fluid 
movement that is influenced by an applied magnetic field. In FHD, 
there is usually no electric current flowing in the fluid, as 
opposed to magnetohydrodynamics (MHD) . {In MHD flows the body force 
acting on the fluid is the Lorentz force that results when electric 
current flows at an angle to the direction of an applied magnetic 
field) . A ferrofluid is a colloidal suspension of solid magnetic 
particles in typically a Newtonian parent liquid such as kerosene. 
In a true ferrofluid, the colloid suspension never settles out 
because the particles are small enough (3-15 nm) that thermal 
agitation keeps them suspended. In addition, particles are coated 
with a dispersant that provides for short range repulsion and 
prevents agglomeration of the particles. A typical ferrofluid 
contains 10” particles per cubic meter and is opaque to visible 
light. Ferrofluids are not found in nature, but are the result of 
laboratory synthesis. Ferrofluids retain their fluid nature even 
in intense magnetic fields (Rosenweig, 1985, and Cowley and 
Rosenweig, 1967). 

The body force in FHD flows is due to a polarization force. 
This requires that the magnetic particles in the ferrofluid align 
in the presence of an applied magnetic field. The particles in a 
colloidal ferrofluid have an embedded magnetic moment. When the 
magnetic field is absent, these particles are randomly oriented, 
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and the fluid has no net magnetization. When ordinary field 
strengths are applied, thermal agitation partially overcomes the 
tendency of the dipole moments to align with the applied field. 
However, as the magnitude of the magnetic field increases, the 
particles become more aligned with the direction of the field. In 
the presence of very intense magnetic fields, the particles may 
align completely, and the magnetization is said to be saturated. 

Applications of ferrofluids are diverse. These include zero 
leakage rotary shaft seals used in computer disk drives, vacuum 
feed throughs, and pressure seals for compressors and blowers. 
Also widely available are liquid cooled speakers that use drops of 
ferrofluid to conduct heat away from the speaker coils. In the 
medical field, the use of a magnetic field can direct the path of 
a drop of ferrofluid in the body, which allows for the directing of 
drugs to a target site. In addition, ferrofluids can be employed 
as a tracer of blood flow in non- invasive circulatory measurements. 
In other areas, high volumes of ferrofluids are utilized in sink- 
float separation processes that use the artificial high specific 
gravity imparted to a pool of ferrofluid subjected to a magnetic 
field. This process is used to separate scrap metals and is also 
used to sort diamonds. Ferrofluids are also being considered as a 
possible candidate for ink in high-speed, silent printers. 

Modeling of physical problems which involve ferrofluids 
typically utilize on of three modeling approaches: 


1) inviscid flow 



2) viscous flow with a symmetric stress tensor 

3) viscous flow with a non-symmetric stress tensor. 

Inviscid flows 

This approach is utilized when there exists no shearing 
stresses in the flow. This is the case when flow configurations 
have no bulk motion or when the flow is irrotational (Cowley and 
Rosenweig, 1967). Such problems include wave motion on a surface 
and/or the study of interfacial processes. 

The governing equations for these flows include: 

The Continuity equation 

V ■ ji = 0 (2.1a) 


The Momentum equation 

+ u • Vu = -A P* + Jl 0 MVH (2.1b) 

where is the velocity vector, u 0 is the magnetic permeability of 
free space, P* is the modified pressure, M is the magnitude of the 
magnetization of the fluid, jj/ and H is the magnitude of the 
magnetic field vector applied to the fluid, a. The modified 
pressure P* is given by 

P* = P( P/ D + \i 0 f* (v^g) BtT dM + \i 0 f*M dH (2,2) 

P* = P(p,D + P a + P a 

The magnetostrictive pressure is given by P a , and the fluid magneto 



pressure by P m . 

If the flow is irrotational, the velocity can be expressed in 
terms of the potential, 4>, with 

li = V<j> (2.3) 

Substitution of this expression for the velocity into the 
continuity equation leads to the Laplace Equation 

V 2 <J>= 0 (2.4) 

The magnetic field is governed by Maxwell's equations, which, 
are the following: 

Vx£--|f (2.5) 

where £ is the electric field and £ is the magnetic induction field 
which is defined as 

i 0 (H + tf) (2.6) 

and 

U - (|i - 1) U (2.7) 

where u is the magnetic permeability of the fluid. 

V • £ = p e (2.8) 

where fi is the electric displacement and p { is the free charge due 
to electricity. 



(2.9) 


V x H = + 


dn 
dc ‘ 


where <i £ is the electric current. 


V- Jf = - 



( 2 . 10 ) 


V • & = o 


( 2 . 11 ) 


If there is no electric field/electric current, then Maxwell's 
equations reduce to: 


v x a = & 


( 2 . 12 ) 


V • = 0 


(2.13) 


It follows that with £ proportional to fl, or for uniform jj, 


V • H = 0 


(2.14) 


Since the magnetic field is irrotational, £ can be expressed 
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as 



a = -V T|r 


(2.15) 


where 7 is the magnetic field potential. Substitution into 
equation 2.14 yields 


V» i|r = 0 


(2.16) 


Typical boundary conditions at an interface and/or free 
surface include: 

The Kinematic Condition 

The kinematic condition relates the interface deflection to 
the motion of the adjacent fluid. When the surface position 
changes with time, the location of the interface can be represented 
as 

z = z 0 (x, y ; t) 

If a Monge function is defined by Fe(x, y, z; t)s z - z 0 (x, y; t) , 
then the kinematic equation at the interface is 

^ + Uj/ V Fe = + u* • V Fe = 0 (2.17) 


Normal Force Balance 

Pi * -|n„ • P2 + -|n 0 + ° V (2.18) 
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where subscripts i and 2 denote fluid properties above and below 
the interface. The surface tension of the interface is denoted by 
a. The outward normal to the interface is given by n = 
VFe/ | VFe| 

Also, other boundary conditions which depend on the 
configuration of the problem are needed. These involve constraints 
on the magnetic field and fluid variables. 

Viscous flows with a symmetric stress tensor 

When considering viscous flows with g and g collinear, the 
viscous stress tensor, 2 , is symmetric. A ferrofluid will always 
have a symmetric stress tensor if it is superparamagnetic, which 
means the solid particles instantaneously align with the magnetic 
field (M and g are collinear) . The equation for the viscous stress 
tensor is given by 


X = - P fi 


H 


du £ 

(, 55 


•ft 


+ a 


a 




( 2 . 19 ) 


The equation for the magnetic stress tensor is 

Tj, = - l/ 0 % 0 (-^f)*,e \ Po W&ij + B i h j (2-20) 


Flow configurations involving incompressible viscous fluids 
and a symmetric viscous stress tensor are governed by the 



continuity equation and the conservation of momentum equation as 
well as Maxwell's equations. In this case, the conservation of 
momentum equation is given by 


P 


DU 

DC 


- vp * pa * n v 2 u - v[|* 0 am 

+ I l 0 lf7H 


( 2 . 21 ) 


The equation assumes that the flow is incompressible and isothermal 
with constant viscosity, i). Here, 0 is the temperature of the fluid 
and u= p' 1 . Also, the body force term of gravity is included in 
equation (2.21) and also could be included in equation (2.1b) if so 
desired. 

Additional boundary conditions arise at the interface/ free 
surface for viscous flow problems. These include the continuity of 
the tangential velocity components at the interface and the 
implementation of a tangential stress balance at the interface. 

Viscous flows with a non-symmetric stress tensor 

Up to this point, the magnetization M has been collinear 
with the magnetic field, &, which is the case in static 
equilibrium. If superparamagnetism in the ferrofluid is obtained, 
collinearity is an adequate assumption, because the direction of m 
rotates freely within the magnetic particle. Conversely, in 
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paramagnetic fluids, the particles in the fluid have a magnetic 
moment that is locked to the orientation of the particle. If h 
changes direction, m responds by particle rotation, which is a much 
slower process because it is resisted by fluid viscous-drag torque. 
The result constitutes a body couple of u-m x h. When this body 
couple is present, the viscous stress tensor is no longer 
symmetric . 

This phenomenon occurs in several cases. For example, if a 
beaker of magnetic fluid is subjected to a rotating magnetic field, 
a non-symmetric stress tensor will result. In addition, any flow 
that possesses vorticity and is subject to a steady magnetic field 
will have some degree of antisymmetric stress present. Also, this 
phenomenon occurs in cases in which a steady field is imposed on a 
ferrofluid moving across a stationary flat plate. The result in 
this third case is that the velocity in the fluid in the boundary 
layer adjacent to the plate becomes coupled to the field. The 
boundary flow is rotational, tending to reorient the magnetic 
particles of the ferrofluid in the magnetic field. Thus relative 
rotation of particles and the carrier liquid results. However, if 
the field is parallel to the plate and perpendicular to the flow, 
there is no coupling, because the particles can rotate freely with 
their axes parallel to the field direction. 

A system may not only exchange linear momentum with its 
surroundings, but it can also exchange angular momentum. The rate 
of linear momentum gained from the surroundings is accounted for by 
the body force g. In addition, the surroundings may transmit 



angular momentum through a body couple, S- A surface traction 
exists 


t,, = fl 'X (2.22) 

per unit area of the surface, (where 2 is the fluid's stress 
tensor), which is either accumulated within the system, or is 
removed to the surroundings as linear momentum. The surroundings 
also exert on the surface of the system a couple su> per unit area. 

Let £ denote the total local density of angular momentum. 
Then i is written 


L ~ Z x U + £ (2.23) 

where s. x u is the external angular momentum and a is the internal 
angular momentum, or spin. The spin field describes the rotation 
of the magnetic particles in the ferrofluid and the viscous fluid 
that is entrained by the particles. 

Coupling exists between internal and external forms of angular 
momentum. The integral form of the equation balancing the total 
angular momentum is 

— f p(£+£xii)c?V=r(p£ + i + pZ)dV + 

DtJv Jv (2.24) 

+ x x CjdS 
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c n accounts for the presence of a surface-couple density acting on 
the surface, and is expressed by the following: 


£n = A • (1 c x + j Cy + £ c z ) = A • q (2.25) 

The expression in parentheses is the couple stress tensor Q. This 
couple can arise from viscous diffusive transport of internal 
momentum. 

Thus, the differential equation governing the change in total 
angular momentum is 


where 


p_£ (2 + z + u) =pfi + £xpi; + v*£ + 
£ x (VD + d 


d - e ijJc T jk 

€ ijk = +1 if ijk = 123, 231, or 312 

-1 if ijk = 321, 132, or 213 

0 if i = j, i = k, or j = k 


(2.26) 
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Chapter III 


Computational Methods For Solving 
Free Surface Flows 

Free surface flows can be divided into three main areas: (1) 

those concerned with wave motion on the surface, (2) bulk fluid 
motion in a configuration which involves a free surface, and (3) 
interfacial processes on the two-dimensional free surface itself. 
Solid boundaries can further complicate the configuration, giving 
rise to finite domain problems involving both free and fixed 
boundaries . 

Computational solutions techniques will change, depending on 
whether the flow is viscous or inviscid, steady or unsteady, and 
whether the analysis is linear or nonlinear. 

The governing equations for free surface/interfacial flow 
problems include: 

The Continuity equation 

+ V • pu = 0 (3.1a) 

or 

V • u - 0 (3 . lb) 


in the case in which the fluid is incompressible. 
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Conservation of Momentum equation 


p(-^+ji , Vji)=-VP+p2+jiV 2 u (3.2) 

for viscous, incompressible flows. If the flow is inviscid, the 
last term on the right side of the equation is absent. 

Boundary Conditions for inviscid flows 

1. If the flow is inviscid, the normal velocity components at an 
interface must be equated. 

2 . The normal force balance for free surface flows is 

AP=oV-J9 (3.3) 

where n is the unit outward normal to the interface, o is the 
surface tension and AP denotes the difference in pressures of the 
fluids above and below the interface. 

P = P 0 ~ p|| - -f V* -V<t> - pgz (3.4) 

on z = 1 + 

where $ is the velocity potential obtained from u = V $, and tj is 
the shape of the interface. 

3 . The Kinematic Condition 

Define Fe = z - n(x,y,t), then 



This ensures that particles on an interface remain on the 
interface . 

Boundary Conditions for Viscous Flows 

1. Velocity components corresponding to upper and lower regions 

must be equated. 

2. The kinematic condition (the same equation as for inviscid 

flows) . 

3. A normal and tangential stress balance must be utilized. 

If the configuration of the problem is temperature dependent, 
then the Conservation of Energy equation must also be utilized. 

pc p || = k V 2 0 (3.6) 

where 8 is the temperature, k is the thermal conductivity of the 
fluid, and c p is the specific heat of the fluid at constant 
pressure. Note that contribution from viscous heating and any 
potential thermal source have been neglected 

These above equations illustrate in general the system which 
must be solved. However, additional thermal and fluid boundary 
conditions specific to the particular configuration would be needed 
in order to formulate the problem properly. Once formulated, the 
problem would be solved via techniques of computational fluid 
mechanics . 

From the literature search, a database was formed by culling 

articles published in Journal of Computational Physics and 

International Journal of Numerical Methods in Fluids and further 
references obtained from selected journals. This database is 
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included in the report as Appendix E. From these references, it 
was determined that the type of computational grid used in free 
surface problems can be either Eulerian or Lagrangian in nature. 
In an Eulerian grid, the cells remain fixed in the domain of the 
problem and fluid movement through the cells is tracked. In 
contrast, Lagrangian grids have cells that move along with the 
fluid so that the same cell in the mesh is associated with the same 
fluid element. The choice of using either method is application 
specific and both configurations are widely used. 

Eulerian Methods 

It has also been discovered that Eulerian methods can be 
further divided into three categories: 

1. Fixed Grid Methods 

2 . Adaptive Grid Methods 

3 . Mapping Methods 

In the fixed grid method, the grid does not change in the 
domain of the problem. Finite difference methods have been mainly 
used for fixed grids because these are predominantly utilized for 
the simplest types of flow configurations. The two basic ways of 
tracking the interface are surface tracking methods (Marker and 
Cell) and volume tracking methods (Volume of Fluid) . These methods 
are detailed in Appendix A and B. 

Surface tracking involves specifying a set of marker particles 
on the surface. The particles move according to the governing 
equations of the flow — usually the Navier-Stokes equations and the 
continuity equation. In volume tracking, however, each interface 
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cell is specified. by a fractional volume of fluid in the cell. A 
function is defined that denotes which cells have fluid, which ones 
are partially filled, and which ones are empty. The surface cells 
are those that are adjacent to empty cells. 

In general, the moving interface does not coincide with a grid 
line in fixed grid methods. In this case, special book-keeping 
procedures must be included in the algorithm in order to handle 
this. It is difficult to calculate the position of the interface 
accurately. 

Some advantages to using fixed grids are: 

1. It is possible for interfaces to undergo large deformations 
with little loss of accuracy. 

2. It is straightforward to handle multiple interfaces. 

Adaptive grid methods involve moving meshes where the motion 

of the grid is linked to the deformation of the interface. They 
are still considered Eulerian because the domain remains fixed. 
These methods are much more versatile in that either steady or 
unsteady flow can be handled. The main reasons for using adaptive 
grids are: 

1. The grid can be moved in the interior of the solution domain 
so that local accuracy can be attained, if desired. 

2. In moving boundary problems, the grid can conform to geometry 
changes in the problem. 

Finite elements are preferred for adaptive grids because they 
are more efficient in handling complex geometrical configurations. 
A disadvantage of this method is that if there are large 
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distortions in the free surface, the grid can become twisted, 
causing the solution to the field variables to be highly difficult. 

Further exposition of adaptive grids can be found in Appendix 

C. 

In the mapping method, an unknown, irregularly shaped flow 
domain is transformed onto a fixed and regularly shaped 
computational domain. By utilizing this technique, it becomes 
fairly easy to obtain a finite difference mesh for which boundaries 
coincide with mesh lines. The mapping function is an explicit 
unknown and must be determined along with the field variables. 
However, it is quite difficult to determine, in the general case, 
what the appropriate mapping function is. Appendix D presents some 
details on mapping methods for the specific problem of a bubble 
moving in an unbounded liquid. 

Lagrangian Methods 

Lagrangian methods are well suited for free surface problems 
because they allow material interfaces to be specifically 
delineated and precisely followed. In addition, interface boundary 
conditions can be easily applied. 

Lagrangian methods can be categorized in two ways: 

1. Strictly Lagrangian methods 

2 . Lagrangian methods with rezoning 

In strictly Lagrangian methods, the computation grid topology 
is fixed and it moves along with the fluid. Clearly, these cases 
should only be used for limited evolution times so that mesh 
tangling does not occur. 
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A well known Lagrangian scheme of type 1 is a method proposed 
by Hirt et al (1970), referred to as LINC (Lagrangian 
Incompressible) . This technique is used for transient two- 
dimensional flow of viscous and incompressible free surface 
configurations. This method employs a finite volume discretization 
with quadrilateral cells. Positions, velocities, and body 
accelerations are defined at the vertices. The pressures and 
stress tensors are stored at the center of each cell. The pressure 
of a given cell is derived by the constraint that the volume of 
each cell remains constant when vertices are moved to a different 
position. As a result, a Poisson-like equation is developed. 

Lagrangian methods with rezoning are introduced when the 
computational mesh become severely distorted. In this case, a new 
mesh is developed and information from the old mesh is transferred 
to the new one. Typically, this procedure calls for a Lagrangian 
phase followed by a rezoning phase where mesh points are moved to 
the prescribed positions. If rezoning occurs at every time step, 
the process is called continuous rezoning. However, if rezoning 
occurs after many time steps, but before the mesh becomes tangled, 
the process is called general rezoning. 

Rezoning of the Lagrangian boundary vertices simplifies 
further the treatment of strongly distorted interfaces. This 
involves tangential rezoning, which equalizes the spacing of 
vertices along an interface, and normal rezoning, which eliminates 
the interface distortion that cannot be resolved with the employed 
grid structure. Newer Lagrangian methods can track an interface 
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through large distortions with the use of rezoning even with 
surface tension effects included. 

A further issue raised in the literature concerns the efficacy 
of finite difference techniques versus finite element methods. It 
is noted that the choice of method depends on the configuration of 
the problem. Finite differences are preferred when fixed grids and 
mapping methods are employed. Finite elements, however, are more 
suited to moving boundary problems where adaptive grids are used. 
Also, finite elements are more efficient when complicated 
geometries of the flow exist. 
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Chapter IV 


Conclusions and Recommendations 

The combination of the topics of ferrohydrodynamics and free 
surface fluid dynamics proves to be highly relevant to commercial 
research applications. There also exists a considerable amount of 
new computational methods for solving free surface fluid mechanics 
problems, which could be utilized in ferrohydrodynamic 
applications . 

The area of ferrohydrodynamics is a relatively new area of 
fluid mechanics. Use of ferrofluids allows for control of the 
fluid configuration by the applied magnetic field. This has both 
current and potential applications in a terrestrial environment as 
well as a microgravity environment. 

It is recommended that future work in this area focus on the 
formulation of specific research problems which are scientifically 
important and whose results would be of interest to industry. 

It is anticipated that future proposals will be written in 
this area. 
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Appendix A 

The Marker and Cell Technique 


This method is used for problems involving time dependent 
motion of viscous, incompressible fluids in a two-dimensional 
configuration. The fluid can be bounded in part by walls of an 
irregular box or by lines of reflective symmetry. Initial and 
boundary conditions are supplied. 

The solution methodology utilizes primarily finite difference 
approximations applied to Naiver-Stokes equations . The dependent 
variables are pressure and two velocity components. 

The finite differences apply to space and time. The region is 
divided into small rectangular cells where field variables are 
given by single, average values. Time changes are represented by 
a sucession of field variables separated by small increments of 
time . 

The complete velocity field in known at the beginning of the 
cycle. The coordinates of the marker particles are known (these 
show which regions have fluid and which do not) . The pressures are 
found such that the rate of change of velocity divergence vanishes 
everywhere. The two components of acceleration are found, 
multiplied by time, incremented per cycle, and the change in 
velocity is added to the old values. The marker particles are 
moved according to the velocity components. Adjustments are made 
for passing marker particles across cell boundaries. Velocity 
modifications are made when a cell, previously full, is emptied, or 
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vice versa. The marker particles do not participate in the 
calculation, but they do show fluid element trajectories and 
relative position of fluids. 

The equations utilized are 
Continuity 


du 

dx 


+ 



and 

Conservation of Momentum 


(A.l) 


du + du 2 
dt dx 


+ 


duv + dP 
dy dx 


= v (V 2 u) 


+ 9 X 


(A. 2a) 


dv 

It 


+ 


duv ■ 
~lx 


dy 2 

dy 


+ = v (V 2 v) 


+ 


(A. 2b) 


where, here, P is the ratio of pressure to constant density. The 
fluid is divided into small regions with local field variables. 
They are given indices i and j which count cell center positions 
along the horizontal and vertical directions. 

Cell boundaries are labeled with 1/2 integer values of 
indices. The rectangular cells have dimensions of fix and 6y. 

Superscripts are used to number the time cycle and fit is the 
time increment per cycle. 
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The finite difference approximations for the continuity and 


momentum equations are 



(u J(i ) 2 -(u i+ i fi ) 2 + (v i+ i i+ i) 

2 2 2 2 

— 

, + - 2 ^ 

+ V ( 2 


fix 2 




+ ~ 2 ^4**, 


fiy 2 

+ p i,J ~ p i*i,j 

fix 


+ S’* 


y i>J y l,j _ 3 _ 2 2 J 2 2 J 2 2 ,J 2 

fit fix 

+ (Vi.i ) 2 - (Vi .^) 2 

fiy 

- 2 V, 


+ v ( 


v i*i.j+\ + ‘ * v i.j4 

fix 2 

^4 + ^4 

fiy 2 
+ p i,i ~ 

fiy 


(A. 3) 


(A. 4) 
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Define 


'i.i 




5x 


<^4 - ^4) 

5y 


(A. 5) 


So the continuity equation is 


D ifJ = 0 


(A. 6) 


From the momentum equations we obtain 


"fiF 


D 11 - D Jul =- 0li - I 1 - 2 

1>j fix 2 

- p i,J*i * ~ 2 p i,J 

fiy 2 

+ v ( + 2 Eli 

5x J 

+ g i. J*i * ~ 2 j ) 

fiy 2 


(A. 7) 


n+1 

We want D i(j to disappear for every cell at the end of a cycle. 
This leads to 

p ** hi 1 " 2 Ihl + ?LteL± ELt±J-LIhJ. = - r > . (A. 8) 

fix 2 fiy 2 i,:/ 
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where 




= 0: 


'i.i 


+ 


D i.j _ v ( D i+I.j * ElzU. 2 Bid 

fit fix 2 

~ 2 P i,J ) 

5y 2 


(A. 9) 


and 


Qi,s = (^W 2 + - 2 + (Vi.j+i ) 2 + Vi.i-i ) 2 - 2 

21 , + i) (v j+i .,i) + (C7,_1 (V, i , i)] 

fix 5y 

2 t (~ U l+±,j+ l) + ( U l-±,j+l) 

+ 2 ,J 2 2 ,J 2 2 ,J 2 2' J 2 

fix 6y 


(A. 10) 

First, R t (j is found for every cell, then Pi^ is obtained for 
the cells. Finally, the new velocity components are found from the 
momentum equation. 

To account for the free surface, the new position of the 
surface must be calculated. This calculation can be done because 
of the coordinate system of marker particles which follow the 
motion of the elements throughout the fluid. The marker particles 
are placed in cells with fluid, and they move with the local 
velocity. Cells with no marker particles are said to have no fluid 
in them and cell with marker particles that are adjacent to empty 
cells are called surface cells. 

Boundary Conditions 

There two types of walls: free slip and no slip walls. For 
a free slip wall, normal velocity reverses and tangential velocity 
remains the same. However, for a no slip wall, the opposite is 
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true. Also, in a free slip wall the pressure is 


p' = P ± g x fix 


(A. 11) 


for a vertical wall, and 

P' = P ± g y by 


(A. 12) 


for a horizontal wall. 

For a no slip wall, 

P' = ± g x bx ± (— (A. 13 ) 


for a vertical wall, and 


P' 


= Px 


± 


sr y ± ( 


2w m 

~w ] 


(A. 14) 


for a horizontal wall. 

At the free surface, D = 0 for surface cells. The pressure 
boundary condition is derived from either the disappearance of the 
normal stress component or from equating it to the applied external 
pressure. 
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Appendix B 

The volume of Fluid (VOF) Method 


This method involves the introduction and use of a function F 
that has a value of one at any point occupied by a fluid and zero 
otherwise. The average value of F in a cell represents the 
fractional volume of the cell occupied by the fluid. Cells with 
values between one and zero must contain a free surface. The 
normal direction of the boundary lies in the direction in which F 
values change most rapidly. 

The time dependence of F is governed by the equations 


dF 

dt 


+ u 


dF 

Sc 


+ v 



(B.l) 


F moves along with the fluid, and it is the partial differential 
equation analog of the marker particles. In this case, F is a flag 
identifying cells containing fluid. 

The method follows a region of fluid rather than surfaces, 
has minimum storage requirements, and there are no problems with 
intersecting surfaces. 

If the flow is incompressible, then the conservation of mass 
equation 

l£ + !? = o (B - 2) 

must be implemented. If limited compressibility is desired, then 
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will hold, where C is the speed of sound in the fluid. 

For specificity, consider the conservation of momentum 
equation to be given by the Naiver-Stokes equations . 

Discrete values of dependent variables, including F (the 
fractional volume of fluid) are used in the VOF technique. The 
free surface cells are the ones which contain a nonzero value of F 
and have empty cells neighboring them. 

There are three steps for advancing a solution through one 
increment of time, fit. First, approximations of the Naiver-Stokes 
equations compute the first guess for new time level velocities 
using initial conditions or previous time-level values for all 
advective, viscous, and pressure accelerations. Secondly, the 
pressures are adjusted in each cell and consequent velocity changes 
are added to velocities computed in the first step. This is done 
to satisfy continuity. Finally, F must be updated to denote new 
fluid configurations. These steps are repeated to advance a 
solution through any desired time interval. At each step, boundary 
conditions are imposed at all mesh and free surface boundaries. 

This method improves on the old MAC method by using a more 
accurate form of the continuity and momentum equations. In MAC, 
continuity and momentum equations were combined so that the 
convective flux term could be written in divergence form (V 3i Si 
instead of ji • V 34 ) . This ensured conservation of momentum in the 
difference approximations. This does not, however, ensure accuracy 



in a variable mesh. It is better to use the form for advective 
flux. 

An iterative process is used to adjust velocities and 
pressures in each cell in order to satisfy continuity. In each 
cell with fluid, the pressure is changed to either push out fluid 
or to bring it into the cell. This must be dome in a number of 
passes through the mesh. In cells containing a free surface, the 
pressure is assumed to be specified. The surface cell pressure is 
set equal to the value found by linear interpolation between the 
pressure at the surface and a pressure in the fluid. 

In the VOF technique, the boundary of the surface cells is 
approximated by a straight line through the cell. The slope of the 
line is determined and it is moved across the cell to a position 
that intersects the known amount of fluid in the cell. 
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Appendix C 

Adaptive Grid Methods 


Adaptive grid methods involve moving meshes where the motion 
of the grid is linked to the movement of the interface. To 
demonstrate this process, consider the problem solved by Cuvelier 
and Driessen (J. Fluid Mech., 1986) involving the thermocapillary 
free boundaries in crystal growth. The governing equations are the 
Naiver-Stokes equations coupled with the heat equation. Boundary 
conditions include the normal and tangential stress conditions, 
continuity of thermal flux, conservation of volume, and the 
condition that the normal component of velocity is zero at the free 
surface. The goal of this problem is to find the shape of the free 
boundary, <|>, a velocity vector, u, and a pressure field that 
satisfies all of the equations describing the configuration. The 
boundary condition that is relaxed is 

{Oh Re) 2 o y = - B 0 <|> (C.l) 

A 

which is the normal force condition. Oh is the Ohnesorge number, 
equal to 

oh = 14 — _ 

(Po Yo 

Re is the Reynolds number, equal to 

Re-^ ill 

P 
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Bo is the Bond number, equal to 


Bo 


9 Po L 2 
Yo 


An auxiliary problem is defined by finding the solution 
{u,p,T} satisfying all of the remaining boundary conditions. The 
relaxed boundary condition is used to adjust the function 4>. 

The function $ will now be compared with the free boundary of 
the static, problem, <J>. . If f = <J> + then f satisfies the 
following equation 


( — r ) / + B 0 W = G 

[1 + <|> 2 .] 1 


(C .2) 


f • ( 0 ) = f ( 1 ) = 0 


with G given by 


G - - {Oh Re) 2 (o v .(p 0 u)) - f 1 a v {p 0 u)dx x 

Jo 




[i + ♦;'] * 


(C.3) 


The algorithm is based on the following: the free boundary $ 
is iteratively approximated by a sequence <fr°, (fr 1 , . . .(fr 1 defined by 

(i) the static free boundary 

(ii) assume , . . . are known 


A33 



(iii) solve the auxiliary problem with (fc 1 ; the solution is denoted 
by {u 1 , p 1 , T 1 } 

(iv) Solve problem (2) with G 1 = G^^u 1 ,p* ,T*) 

The solution is denoted by Jr 1 . 

(v) <J) 1+1 = 

The auxiliary problem is nonlinear, so it must be solved 
iteratively. The Newton Raphson method is used to linearize the 
nonlinear terms. 
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Appendix D 

Mapping Methods 

Mapping methods make it possible to transform an irregularly 
shaped flow domain into a fixed, regularly shaped computational 
domain. This makes it easy to implement a finite difference mesh 
for which the boundaries coincide with a mesh line. This method 
also maintains a sharp resolution of the interface. 

The general algorithm for mapping methods is as follows : 

1. The governing equations of the flow are introduced and 
boundary conditions are implemented. 

2. An equation for the surface is set up. This is regarded as 
an unknown function. 

3. The equations are non-dimensionalized. 

4 . The coordinate transformation occurs . This solves the 
problem of the difficulty of having an implicit unknown function 
of the surface. It also insures that the boundary lies on a 
coordinate line. The governing equations and the boundaries are 
solved under the assumption that the surface is now a known 
function. Next, the pressure is restored on the basis of the 
calculated stream function of velocity function. When all of the 
other functions are known, the calculation of the original free 
surface is the final step. 

An example of this method is given by Christov and Volkov 
(J. Fluid Mech., 1985) in their study of the steady viscous flow 
past a stationary, deformable bubble. 
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Consider an unbounded volume of viscous incompressible liquid 
having a velocity U. at infinity. There is a resting bubble in 
the flow. 


The 2-D Naiver-Stokes equations in spherical coordinates are 
introduced in a stream function-vorticity formulation. 

The first boundary condition states that at infinity, liquid 
moves with a uniform velocity U. in the z direction. 

Y(r, 0)--| U m r 2 sin 2 0 (D.la) 

4b 


5-0 at r-«> 


(D. lb) 


The other boundary conditions are 


Y(r, 0) 


- -i U m r ? sin 2 0 - f (1 + cos©) x 
2 4nV m 

[1 - exp(-^(l - cos0))] 

P 


(D . 2 ) 


5(r, 0)- - 


&VU. sinO 
4 nv 2 r 


(1 + -~-)exp[i^(l - cos0) ] 


(D. 3 ) 


g = acceleration due to gravity, V = volume of the bubble 


r = *(0) 


Let the equation of the surface be 
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(D. 4) 



Next, the equations are non-dimensional ized. Because r = 
R(0) is unknown, many obstacles arise in solving the problem. To 
overcome this problem, the coordinate transformation takes place. 
The transformation used to simplify the problem is 

x = ti i?(0) , 0=0 (D . 5 ) 

These are substituted into all the equations. A Finite 
difference method can now be used to solve for the function R(0) . 
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Appendix E 


Database of Articles Collected 
Pertaining to Free Surface Flow 


Free Surface Flow During the Filling of a Cylinder, 

Abdullah, Z. and Salcudean, M. International Journal for Numerical 
Methods in Fluids, vol. 11, p. 151-168, 1990. 

Key words: 

Free surface flows Filling of molds Flow Visualization 
Mathematical modeling Control volume method 


A High Resolution Godunov-Type Scheme in Finite Volumes for the 2D 
Shallow-Water Equations, Alcrudo, Francisco and 
Garcia-Navarro, Pilar. International Journal for Numerical Methods' 
in Fluids, vol. 16, p. 489-505, 1993. 

Key words : 

Two-dimensional modeling Finite volumes MUSCL approach 
Upwind differencing 


Numerical Solution of the Generalized Serre Equations with the 
Maccormiack Finite Difference Scheme, Antunes Do Carmo, J.S, Seabra 
Santos, F.J. and Almeida, A.B. International Journal for Numerical 
Methods in Fluids, vol. 16, p.725-738, 1993. 

Key words: 

Serre equations Maccormack's Method Solitary Waves 
Sudden releases 


Flair: Flux Line-Segment Model for Advection and Interface 

Reconstruction, Ashgurz, N. , and Poo, J.Y. Journal of 
Computational Physics, vol. 93, p. 449-468, 1991. 

Key words: 

Cell volume fraction approach Line segments Capillary forces 


Adaptive Zoning for Singular Problems in Two Dimensions 
Brackbill, J.U. and Saltzman, J.S. Journal of Computational 
Physics, vol. 46, p. 342-268, 1982. 

Winslow's method Adaptive mesh 


A Continuum Method for Modeling Surface Tension, Brackbill, J.U., 
Kothe, D.B., and Zemach, C. Journal of Computational Physics, vol. 
100, p. 335-354, 1992. 

Key words: 

Surface tension Fluid interfaces Continuum method 
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Modelling the Initial Motion of Large Cylindrical and Spherical 
Bubbles, Bugg, J.D. and Rowe, R.D. International Journal for 
Numerical Methods in Fluids, vol. 13, p. 109-129, 1991. 

Key words : 

Bubbles Finite difference Initial motion 


Application of a Vortex Method to Free Surface Flows 

Chen, Liyong and Vorus, William. International Journal of Numerical 

Methods in Fluids, vol. 14, p. 1289-1310, 1992. 

Key words: 

Vortex method Free surface flows Body-wave interactions 


Numerical Investigation of the Steady Viscous Flow Past a 
Stationary Deformable Bubble, Christov, C.I. and Volkov, P.K. 
Journal of Fluid Mechanics, vol. 138, p. 341-364, 1983. 

Key words : 

Moving Boundaries Coordinate transformation Steady Viscous 
Deformable bubble 


Explicit Streamline Method for Steady Flows of Non-Newtonian 
Matter: History Dependence and Free Surfaces, Chung, S.G. and 

Kuwahara, K. Journal of Computational Physics, vol. 104, p. 
444-450, 1993. 

Key words: 

Finite difference method Non-Newtonian matter History dependence 
Free surface 


An Orthogonal Mapping Technique for the Computation of a Viscous 
Free-Surface Flow, Cliffe, K.A., Tavener, S.J., and Wheeler, A. A. 
International Journal for Numerical Methods in Fluids, vol. 15, p. 
1243-1258, 1992. 

Key words: 

Free surface Viscous flow Finite element method 
Orthogonal mapping 


Thermopapillary Free Boundaries in Crystal Growth, Cuvelier, C., 
and Driessen, J.M. Journal of Fluid Mechanics, vol. 169, p. 1-26, 
1986. 

Key words: 

Two-dimensional flow Thermocapillary flow Steady Crystal 
growth Finite element method 
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A Newton's Method Scheme for Solving Free-Surface Flow Problems, 
Dandy, David S. and Leal, Gary L. International Journal for 
Numerical Methods in Fluids, vol. 9, p. 1469-1486, 1989. 

Key words : 

Bubble Intermediate Reynolds number Newton's method 


An Efficient Surface-Integral Algorithm Applied to Unsteady Gravity 
Waves, Dold, J.W. Journal of Computational Physics, vol. 103, p. 
90-115, 1992. 

Key words : 

Unsteady motion Two-dimensional flow Cauchy's integral theorem 
Numerical instabilities 


Three-Dimensional Streamlined Finite Elements: Design of Extrusion 
Dies, Ellwood, Kevin, R, Papanastasiou, T.C. and Wilkes, J.O. 
International Journal for Numerical Methods in Fluids, vol. 14, p. 
13-24, 1992. 

Key words : 

Streamlined Finite elements Extrusion Free surface flows 
Die design Three-dimensional flows Three-dimesional finite 
elements 


Consistent vs. Reduced Integration Penalty Method for 
Incompressible Media Using Several Old and New Elements, 

Engleman, M.S. and Sani, R.L. International Journal for Numerical 
Methods in Fluids, vol. 2, p. 25-42, 1982. 

Key words: 

Penalty Method Incompressible flow Reduced quadrature 
Finite elements 


Finite Element Method for Time -Dependent Incompressible Free 
Surface Flow, Frederiksen, C.S. and A.M. Watts. Journal of 
Computational Physics, vol. 39, p. 292-304, 1981. 

Key words: 

Finite element method Time dependent flow Vertically moving 
plate Circulation in rectangular channel 


Flux Difference Splitting for Open-Channel Flows, 

Glaister, P. International Journal for Numerical Methods in Fluids, 
vol. 16, p. 629-654, 1993. 

Key words: 

Shallow-water equations Subcritical and supercritical flows 
Open channels 
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A Finite Difference Technique for Solving the Newtonian Jet Swell 
Problem, Han, Chin-Tsu, Tsai, Chin-Chin, Yu, Tsai-An, and Liu, 
Ta-Jo. International Journal for Numerical Methods in Fluids, vol . 
15, p. 773-789, 1992. 

Key words : 

Finite difference method Newtonian jet swell 


Lagrangian Finite Element Method for Free Surface Navier-Stokes 
Flow Using Fractional Step Method, Hayashi, Masahiro, Hatanaka, 
Katsumori, and Kawahara, Mutsuto. International Journal for 
Numerical Methods in Fluids, vol. 13, p. 805-840, 1991. 

Key words : 

Finite element method Lagrangian description Fractional Step 
Navier-Stokes equations Linear interpolation Free surface 


Numerically Induced Phase Shift in the KdV Soliton 

Herman, R.L., and Knickerbocker, C.J. Journal of Computational 

Physics, vol. 104, p. 50-55, 1993. 

Key words : 

KdV equation Finite difference Phase shift Perturbation theory 


Vol.ume of Fluid (VOF) Method for the Dynamics of Free Boundaries, 
Hirt, C.W. and Nichols, B.D. Journal of Computational Physics, vol. 
39, p. 201-225, 1981. 

Key words: 

Fractional volume of fluid Incompressible Finite difference 


A Lagrangian Method for Calculating the Dynamics of an 
Incompressible Fluid with Free Surface, Hirt, C.W., Cook, J.L., and 
Butler, T.D. Journal of Computational Physics, vol. 5, p. 103-124, 
1970. 

Key words: 

Transient flow Viscous Incompressible Lagrangian coordinates 


The Effect of Surface Contamination on Thermocapillary Flow in a 
Two-Dimensional Slot, Homsy, G.M. and Meiburg, E. Journal of Fluid 
Mechanics, vol. 139, p. 443-459, 1984. 

Key words: 

Insoluble surfactants Steady thermocapillary flow 
Interfacial tension Lubrication theory 
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The Numerical Solution of the Thin Film Flow Surrounding A 
Horizontal Cylinder Resulting from a Vertical Cylindrical Jet, 
Hunt, Ronald. International Journal for Numerical Methods in 
Fluids, vol . 14, p. 539-556, 1992. 

Key words : 

Thin film flows Free boundary problems Parabolic equations 


Numerical Methods for Tracking Interfaces, Hyman, James M. Physica 
12D p. 396-407, 1984. 

Key words : 

Front tracking Interface tracking Surface tracking Volume 
tracking 


Numerical Solution of the Eulerian Equations of Compressible Flow 
by a Finite Element Method Which Follows the Free boundary and the 
Interfaces, Jamet, P. and Bonnerot, R. Journal of Computational 
Physics, vol. 18, p. 21-45, 1975. 

Key words : 

Eulerian equations Compressible flow Finite elements 


Finite Element Method for Shallow Water Equation Including Open 
Boundary Condition, Kodama, Toshio, Kawasaki, Tomoyuki, and 
Kawahara, Mutsuto. International Journal for Numerical Methods in 
Fluids, vol. 13, p. 939-953, 1991. 

Key words: 

Finite element method Shallow water equations Open boundary 
Parameter identification Tokyo Bay 


A Total Linearization Method for Solving Viscous Free Boundary 
Flow Problems By The Finite Element Method, Kruyt, N.P., Cuvelier, 
C., Segal, A., and Van Der Zanden, J. International Journal for 
Numerical Methods in Fluids, vol. 8, p. 351-363, 1988 
Key words: 

Navier-Stokes equations Finite element method Viscous flow 
Free boundary flow 


Turbulent Free Surface Flow Simulation Using a Multilayer Model, 
Lai, C.J. and Yen, C.W. International Journal for Numerical Methods 
in Fluids, vol. 16, p. 1007-1025, 1993 
Key words: 

Multilayer model Free surface Recirculating flow 
Curvilinear co-ordinates Non-staggered grid Depth correction 
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Variational Formulation of Three-Dimensional Viscous Free-Surface 
Flows: Natural Imposition of Surface Tension Boundary Conditions, 

Lee-Wing Ho and Anthony T. Patera. International Journal for 
Numerical Methods in Fluids, 1991. 

Key words : 

Curvature Finite Element method Free surface flow 
Navier-Stokes equations Spectral element method Surface tension 
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SECTION B 



Interface Response of the Finite Length Liquid Column 

to Longitudinal Forcing 
and 


Interface Stability of the Infinite Liquid Column 
to Forcing Normal to the Column Axis 



Axial forcing of an inviscid finite length fluid cylinder 

M. J. Lyell 

Departmem of Mechanical and Aerospace Engineering, West Virginia University. 
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(Received 8 August 1990; accepted 7 March 1991) 

Current interest in microgravity materials processing has focused attention upon the finite 
fluid column. This configuration is used in the modeling of float zones. In this Brief 
Communication, the incompressible inviscid finite length fluid column is subjected to an axial, 
time-dependent disturbance. The response of the fluid system and the resulting interface 
location are determined via Laplace transform methods. 


The presence of vibrations will affect the fluid dynam- 
ical environment in which materials processing occurs. Ex- 
periments performed on Spacelab D1 investigated the ef- 
fect ofg-jitter on the stability of finite fluid columns. It was 
found that long columns are particularly sensitive to resid- 
ual accelerations and impulses. 1 The analytical effort has 
been concerned with the axial forcing of slender liquid 
bridges. 2 ' 5 In these analyses, a one-dimensional “slice" 
model has been used. Such a restriction is invalid in the 
situation in which the liquid column is not slender, yet 
whose height and radius are within the static stability lim- 
its for fluid columns. The work of Sanz 6,1 has addressed the 
problem of the unforced modes of oscillation of a nonslen- 
der finite fluid column. 

In this work, the problem of the axially forced finite 
fluid column is investigated, with the axial forcing consid- 
ered periodic in time. The analysis is linear, inviscid, in- 
compressible, and axisymmetric. It utilizes Laplace trans- 
form techniques. The direction of the imposed body force 
is taken parallel to the longitudinal axis, and in the nega- 
tive e z direction. The undisturbed column is cylindrical, 
and is surrounded by an inert gas. 

Quantities will be nondimensionalized as follows: 


panded in terms of Fr about a state of zero mean motion 
This results in the linearized system of governing equations 

V-u = 0, i 3a ) 


d(V6) 

— V— = - Vp -gU)e r (3b) 


with u = Vd. For definiteness, g(t ) is chosen to be sin(f). 

The boundary conditions at the upper and lower disks 
are that the normal component of the velocity be zero. 
Also, an anchored triple-contact line assumption at both 
disks will further restrict the interface behavior. These are 


d<(> 


dz 


7=5 A 


= 0 , 



(4a) 


/(A,r)=/( - A,r)=0. (4b) 

The constant A acts as a slenderness parameter, and is 
given by (L/2R). As the fluid is incompressible, volume 
must be conserved, and so 



z,t)dz= 0. 


(4c) 


RI = x, cof r 7=r, Rco/G=u, pajjR*p=p, (1) 

and R is the radius of the column. The forcing frequency is 
a>f. Pressure and velocity fields are indicated by p and u, 
respectively. The density of the column is given by p. A 
tilde indicates dimensionless quantities. 

Equations are to be rendered nondimensional through 
the use of ( 1 ). This results in 


Because of axisymmetry 




(4d) 


At the fluid colutnn/inert gas interface, both the kine- 
matic and normal force balance conditions must be im- 
posed. These are given by 


V-u = 0, (2a) 

du 

j { + u-Vu = -Vp- Ft g(t)e p (2b) 

and the nondimensional interface equation is 

f t — r ~ 1 —Frf(z,t) (2c) 

with f(z,t) representing the functional form of the pertur- 
bation to the interface. The tildes have been dropped for 
convenience. The nondimensional parameter Fr is a 
Froude type number, and is given by ( G^Ru >}). It is taken 
to be a small parameter. Pressure and velocity are ex- 


-df 34> A 
dt + dz~ 0. 

(4e) 


(40 


and are evaluated on r= 1 (A p indicates the change in 
pressure across the interface). The nondimensional param- 
eter B 0 is given by (,pR } w}/y). Gamma denotes the surface 
tension. Note that Eqs. (4e) and (4f) have been linearized 
[and are at O(Fr)]. 

Both the system of equations ( 3 ) and the boundary/ 
interface conditions (4) are Laplace-transformed with re- 
spect to time. The time dependence is replaced by a pa- 
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ramctcr j. The initial velocity perturbation is taken to be 
zero. The initial position of the interface is taken to be that 
of the undisturbed cylindrical column f 0 =» 1. This is real- 
istic for the case of a zero mean gravity, which is the case 
in the present investigation. This yields 

,u= -▼/*-( l/^ + nft. (5®) 



(5b) 


The uppercase representation indicates the Laplace trans- 
formed variables. 

The velocity potential is found to be 


<t»(w)= 2 <*.U)/o(V)«*(M* -r A)]. (6) 

«-l 

with k m s (ntr/2A). Here. indicates the unknown 
coefficients and /© represents the modified Bessel function. 
Use of the kinematic condition yields 


m 

2 ^.(i)Mi(*,)cos(*,(* + A)) on r=l. (7i) 

..I 

However, use of the normal force balance condition yields 
a forced ordinary differential equation for F(zj) at r= 1, 


+ 2 4>t4.(j)'o(*.)ccsl*.U + A)l 

Mm | 

(7b) 

Note that s appears as a parameter. Also, K(s) is a eon* 
slant of integration. The solution of Eq. (7b) is 

F(zj) »«(*) cos(z) + 6(j)sin(z) 

+ Bolds) + ((777)) * 

£ / i ^„( r )/ o(*.)«»(*«U + A ) J \ 

+ .?,M — — )• 

(7c) 

Both forms given for F(xj) must be consistent, and 
must allow for the remaining boundary conditions to be 
satisfied. In order to determine the forms of the coefficients 
“r(j)" and “b{s),~ it is advantageous to expand the func- 
tions sin(z), cos(z), and (*) in term of cos(k,(z 
+ A )]. These are found in the literature; 0 

^=[2sin(A)/A]{j/l(l - kiJkiJtkzm) 

- VoUl . J ' MJoU ). ( 8 ) 


I 


with 




A 

l j| 


(f>:. 




(9) 


Dim . I = ( ( 1 — |/|(*I* . |) 


— ff»/o (*!». |)s*l- 

The conservation of volume requires that 

a(*)[sin(A)/A] + BoKis) = 0 . ( I 0 t) 

Finally, the anchored triple-contact line condition is 
imposed. This yields a forced set of algebraic equations in 
the coefficients "a(*r and of the form 

(oil )o(x) + ia\2)b{s) ® cl. ( 10b) 

(o21)o(i) + (o22)6(j) =c2, (10c) 

with ol 1 =o21, o22= — a 12, and c2» — cl. and with 


/nt /At *W A *) ,/ 2 sin(A)\ 

(flll)= cos(A) + ( — - — J 

V Vo C*ii) 

* Kl-tUlOt-l' 

M 2 >- -«.<*> 

r Vo (^ i «* i ) 

-.0 l ( l - kj m .,) Dim+,r 

y fyo(kjm + t ) 1 

-.0 (*L.A— 1>| * 

and with 

*>,«[*,(» -*})/,( V - Vo(V^]- 

It is found that 


(10d) 


(10e) 


(100 


a(s) =0 and 6(s) «■ (d/nl2). (lOg) 

Recall that the expression for Au, contained the factor 
a(s). Thus A}* will be identically zero, and only the odd 
modes will contribute to the motion. This reflects the na- 
ture of the forcing, i.e., "f is in odd function with respect 
to the mean plane between the disks. 

The form of the Laplace-transformed perturbation in- 
terface F(zj) is given by 

/(*.*) =6(j)sin(z) + Bo ( ( 777 )) * 


y *Aim.* \is)Bolo{k lm , i)cos[lt, w . ,(z + A)J 

«-o d - *!■♦■> * 

(ID 

with the Ajm . , as previously determined. Taking the in- 
verse Laplace transform of Eq. (11) will yield the expres- 
sion for f(zj). 

In order to gain insight into the form of the perturba- 
tion interface in the presence of a periodic body force di- 
rected along the longitudinal axis of the column, the series 
will be severely truncated. Only one term in each series will 
be retained. The inversion will then be done analytically. 



f:i r, 



T 0.2S 1.25 2.25 3.25 


FIG. I. Deformation f{z,i) versus time, for lime = 0.25, 1.25, 2.25, 3.25. 
£ 0 = 0.50 and A =0.75. 



Then 


/(z,0=sin(r)h(f) + R 0 (z)$' n (/) 4 -cos(* 1 (r + A)] 
X (xTT^T)( a Hr) S inh[a(,-r)]* 


') 


+ b(t) + cos[*,U + A)] 


S) 


x(J sin(r)sinh[a(/ — t))<1t + a sin(f)j . 

° ( 12 ) 


with 


a 


2 


s 2 = 



b(t) = Jf~'(b(s)), 

and (3 2 8 2 = a 1 . 


The inverse Laplace transform for F(zj) [see Eq. 
(11)] was taken numerically for different values of 
Graphical results are presented in Fig. 1. for the case 
(2?„ = 0.5, A = 0.75) at successive (nondimensional) 
times of 0.25, 1.25, 2.25, and 3.25. The maximum value of 
f(z,t) at each of these times is 0.0019, 0.0145, 0.0131, and 
0.002, respectively. The value of the parameter A indicates 
that the column is not slender. Note the sinusoidal form of 
the interface deformation. This is observed for the case of 
the axially directed forcing. 1 

The equation for F(zj) [Eq. (11)] can be viewed as a 
transfer function. Setting s = ico, values of F may be plotted 
in frequency space. This is done, and the results are shown 
in Fig. 2 for A = 2.60, B 0 = 1.00, and A = 1.25 with 
B 0 = 1.00 and 2.00. The curve for the slender column 
(A = 2.60) may be compared with a similar type graph 
found in Ref. 4 (Fig. 3 in that work). It is seen that a spike 
in both graphs occur about the resonant frequency at 
(o = 0.34. The long column is shown to be sensitive to 
lower-frequency disturbances (<u< 1), as indicated by the 


FIG. 2. Transfer function F(z*,u), e* = - A/2, for the cases (A = 2.60. 
8,= 1.00), ( A = 1.25. B 0 = 1.00). and (A = 125. B 0 = 2 .00). 


relatively large functional values [for F(z,s)] in this work 
and in Ref. 4. Note that the numerical value of F(z*,a>) 
( with z* = — A/2 ) compares well with that of the defined 
transfer function perturbation amplitude in Ref. 4. 

It was found that nonslender columns are not similarly 
sensitive to these low-frequency disturbances. This is 
shown graphically in Fig. 2 for A = 1.25, B 0 = 1.00, and 
B 0 = 2.00. Note the low numerical values of F(z*, o) as 
compared with that of the case A = 2.60, B 0 = 1.00. 

The work of Sanz 6 determined the nondimensional 
natural frequency <5 u , of the fluid column for a range of A 
values. In that work, time was nondimensionalized by 
{pR}/Y) xn . It can be shown that this results in the follow- 
ing equivalence (< o}/a >*,,) = R 0 (l/S^,). 

The perturbation interface was determined (numeri- 
cally) in a number of cases for which the (2?o,A) values 
corresponded to eigenvalues of the unforced system. In 
these cases, the inviscid finite length fluid column was 
found to be very sensitive to forcing. 

Various orders of magnitude for the vibration level 
have appeared in the literature. 8 The range of interest of 
parameter B 0 can be constructed for fluid cylinders of a 
given density and radius. The interface shape can be ob- 
tained for specific values of B 0 and A. 
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Introduction 

T HE float tone configuration is used in crystal growth. It 
may be modeled as a liquid column held by surface ten* 
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Sion forces between imp end dTiki In the case of crystal 
growth, the thermal and solutal fields as well as those of 
velocity and pressure are needed to characterize the physics. A 
compelling reason for crystal growth experiments in a micro- 
gravity environment such as that onboard the Space Shuttle is 
that buoyancy forces are greatly reduced . 1 However, this envi- 
ronment is not quiescent due to the presence of impulse type 
disturbances from small thruster firings as well as those from 
periodic vibrations. Given certain “environmental" condi- 
tions. such as the presence of a periodic acceleration field, it is 
possible that a fluid dynamical instability would develop. This 
would adversely effect the crystal growth. 

The crystal growth environment cannot be separated from 
the fluid dynamics of the liquid column. This latter topic is 
the focus of the present work. The question of column inter- 
face stability in the presence of a periodic acceleration field 
having a component normal to the longitudinal axis of the 
isothermal cylinder is investigated. The fluid co!umn is taken 
to be infinite in length. Floquet theory is used in the stability 
investigation. 

Previous work has determined the natural oscillations of the 
liquid column. This work was done first for the case of the 
infinite column 2 and later for the finite length case. In the 
Utter work, both axisymmetric and nonaxisymmetric oscilla- 
tions were considered . 1 * 4 Results for the infinite length case 
were found to be good approximations to those for the finite 
length column, both numerically and with regard to trends. 

Interface behavior of the finite length liquid column in the 
presence of time-dependent forcing has been investigated for 
the case in which the forcing was parallel to the longitudinal 
axis of the column. w The investigations considered the 
column behavior subject to a unit) forcing for both the in vis- 
cid case of general atptti ratio (within static stability Saits 2 ) 
and the viscous case m the slender cotuaa Emit. 4 

The problem of interface a ability of the fluid column in the 
presence of a periodic acceleration field that has a component 
„ normal to the longitudinal axis of the column has not been 
* investigated. Previous work has considered the interface sta- 
bility of a highly idealised infinite stab-like configuration in 
the praenoe of a periodic ncocicrntioo field oriented normal to 
the interface/ Interestingly, this study was motivated by ex- 
periments in mkrogravity. 

Use of the infinite length configuration in this stability 
study results in a simplification in that the standard boundary 
conditions at the sobd end disks are not appfied, and the focus 
lemairu on the interface stability. If an extension of this work 
to the finite length configuration is of interest. Floquet analy- 
sis would be appropriate, although the implementation would 


a****. 


Forauletioo 

The bode eonflguratioo i* that of «a iaflafeo fluid column of 
circular cron section. Tha fluid is iacompmabi*. and tha 
mnouadini medium it of acgflgibk density. Perturbations 
are taken to be brotatiooaL Tb* analysis is Bacar sad inviacid. 
wish the w o n d imtnt io nt B wd towainq equations those of 
cos uia e i ty and conscnrtiioo of flsosncwaa (notarised Eukr). 

Tha frequency of the periodic foronj b denoted by uy. 
Pressure and uetodty fields ire given by p and m. nondiraen- 
sionaflzed as foOowt: 

W*r */'/• r flw/d«a 0W *r?A m p (l) 

lUdcs indicate norid* mf*ri 5 ?wi quantities. 

The continuity and Euler equations are then 

(2a) 

-ifi-fr cos(r)v[(l - S)ria *1 (2b) 

with A> • (Go/ItJf) a Frauds type amber. Go is the ampli- 
tude of the periodic acceleration (kid. The functional form of 
the time-dependent Cortina is selected to be cos ii). 
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The mean state is one of tero velocity; however, the mean 
pressure is time dependent. Consider a wave like perturbation 
propagating on the interface. The governing equations for the 
perturbation are developed as follows. Expand in the small 
perturbation parameter t 

fi = Pmittft) + *P * = <11 = « V* ( J) 

Substitution of Eq. (3) into Eqs. (2) yields the mean system 

= -/>cos(f)Vl(l - /-)sin01 (4) 

(with the parameter Fr of order one) and the (order <) pertur- 
bation equations 

V 2 $ = 0 (5a) 


a(VO) 

3/ 


= - Vp 


(5b) 


The spatial dependence of the velocity potential can be deter- 
mined via solution of Eq. (5a)* which is 

$(f. a, Z , t) = ZA(t)I m (kr)cxp(ikz)txp(im$) (6) 

J m is the mth modified Bessel function. Clearly, the solution 
involves a superposition of the azimuthal modes. 

It is through the boundary conditions that the free interface 
can be determined and the stability characteristics investi- 
gated. Let the equilibrium interface be given by 

Ft o r - 1 - < 9 ( 0 , z» 0 * r - 1 - *ZC(t)txp(ikz + im$) (7) 

Use kinematic condition and the normal force balance at the 
free interface must be satisfied. In addition, the requirement 
of conservation of mass, which reduces to a’ conservation of 
volume condition, must hold. 

The linearized kinematic condition (at order <) is 


(Sa) 


A 

ii r > I with u, = d 4 >/dr. This results in 

E^j^expO** + imtf) » £|/.(*)l'A(/)exp(/*z + imtf) (8b) 

The normal force balance requires that the difference in 
pressure across the interface be equal to the curvature multi- 
plied by the surface tension force. In nondimensional form, 
this is given by 

Bo x A( p —.' + «p) ■ 7 • a (9*) 

with • - vft/lvFrl. Bo is the nondimensional parameter 
{pR'vj/y), with y the surface tension and p the density. At 
order «, the linearized form of this interfacial condition can be 
expressed as 


1 + «ttt + <ht + [Bo co $(/) sin = Bo 


/«u\ 

\9t) 


(9b) 


Fr. required to be of order one, has been set equal to unity- 
The subscripts indicate partial differentiation. This yields 

E((l - m 1 - k*) + Bo cos (r) sin d| x C(r)exp(ikz + 

= LA (r)/,(*)exp(i*z + imt) ^ 

The sin d dependence can be re-expressed as an exponential 
function. Then 

C.(/KI - m 1 - **) ♦ (y)(-0 cos(/)lC.. ,(/) 


$n J a.\*r mw'Or L >c of tqi a j 

(9d) a nonauior.omous scconJ or Jer equation in C.lO, 

with mode coupling (as indicated by the subsenpu m) 

Cm - Id - nt l - k l )/Bo\ \ /«//« )C* 

“ |/;/2/ - .lcos(/K-i)lC.. l (0-C.* l (/)I (10) 

Keep in mind that cos(/) can be rewritten in exponential form. 

It is at this juncture that Floquet analysis is applied. For 
convenience, let E m {() = |dC„(/)/df). Then take 

(C„(f), f.(r)| = E,;: JC m .„ E mJ 1 X «xp{(X + t/)/| (II) 

The constant coefficients C mJ and E mJ are unknown. The 
Ftoquet exponent is denoted by the eigenvalue X, which is in 
genera] a complex number, and which is also unknown. The 
nonautonomous differential system is thus transformed into a 
homogeneous algebraic system for {C mJ , E mi ) with the un- 
known parameter (eigenvalue) X. If Real(X) is greater than 
zero, the interface of the cylindrical column is unstable to the 
growing wave-like perturbation. Of course, the milieu in 
which this disturbance is propagating includes the periodic 
base state pressure. 

Use of Eq. (11) in Eq. (10) (rewritten as two first-order 
modes) yields the infinite algebraic system given by 

(X + U)C m ,i = E mJ (12a) 

(X ♦ = ((I - - k l )/(Bo)\{/m/\m)C mJ 

+ (*/«/* /1»X - i){C m _ j./. i + C*. !.#♦ i 

■* C» ♦ I.#- 1 — C m ♦ i,/ ♦ i) (12b) 

Note that the harmonic modes (indicated by /) as weD as the 
azimuthal modes (indicated by m) are coupled to both their 
preceding and successive modes. 

Several remarks are in order concerning the truncation. 
Once the truncation in m is done, the number of azimuthal 
modes that contribute are fixed. It is to this truncated system 
that Roquet analysis is actually applied. To obtain numerical 
values for X, it is necessary to truncate the number of harmonic 
modes in time, i.c., the range of / values. The eigenvalue 
problem is, therefore, a problem of the truncated system. 

Results 

The results pertain to the eigenvalue solutions of system 
(12a) and (12b). NAG library routines were used in determin- 
ing the eigenvalues. Truncation values of L » 1151, that is, 
- 15 s i & 15, and M ■ 14 were found to be sufficient. Wave 
number values ranged from k *0.10 to k » 3.00. The parame- 
ter Bo was varied from 0.01 to 10.00. 

For * < 1.0, the interface is unstable to the wave-like pertur- 
bation in the presence of a mean periodic acceleration field 
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I, =0.01 


I • • « e « • t • » 
■ 

r • • ♦ • * eeaeeeeee o # • •««♦*•••«« ♦ ♦ 1 
• ♦ • • t • ♦ • e • •* * «t •«««««•«• 1 

.VAV.V.S^.V.V.V.V/AV.V.V.V/ | 

1 

>>: 




1 .%v.\v.;.y. 

1 



4- 


4- 


4- 


4- 


4- 


4- 


10 0.50 1.00 |.S0 2.00 2.50 3.00 1 


Fig. 1 Stability diagram: (be stability of tbe coo fig u ratios to wave- 
Hkc perturbations for a range of Bo parameter values vs wave lumber 
k la shown. The cross-batched area Indicates unstable regions la which 
disturbances are growing In time. The remaining area corresponds to 
that of the marginal stabDtty state. 
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over the range of Bo values considered. (Note that this range 
encompasses four orders of magnitude, see Fig. I.) 

As k is increased, such that 1.00 s k s 1.20, the interface 
remains unstable to the perturbation for the two larger Bo 
values, equal to 10.00 and 1.00. However, marginal stability 
(Real(X) * 0) ensues at the smaller Bo values. A decrease in Bo 
can be interpreted physically as an increase in the surface 
tension. So, as the surface tension increases, the restoring 
force is sufficient to result in marginal stability over this range 
of wave numbers. 

Perturbations corresponding to the larger wave numbers 
(smaller wavelengths) that were considered, 2.0 s): s3.0, 
were found not to grow in time for Bo - 0.01-1.00. That is, 
instability (with Rcal(X) > 0} of the interface to perturbations 
of these larger wave numbers (and smaller wavelengths) occurs 
only for Bo = 10.00; otherwise, Real(X) = 0. An alternative 
physical interpretation to that involving a variation in surface 
tension for differing Bo values can be developed. Since Fr was 
taken to be unity, (the forcing frequency) is proportional to 
Co, the amplitude of the periodic acceleration Held. Utilizing 
this relation in the definition for Bo yields Boa(Go/y). For 
fixed surface tension (and, of course, density and column 
radius) values, an increase in Bo would result from an increase 
in forcing amplitude. It is at the highest such amplitude con- 
sidered that the interface was found to be unstable (with 
Real(X)>0] to the perturbation. 

It is noted that the range of Bo values used corresponds to 
values of Go and w/that would be of interest in a microgravity 
environment for certain ranges of surface tension values. 
(Roughly, lO^j^nnSGo s lO* 2 *,,**, and 0.5 Hz<«/<5 Hz 
for y values of MOO dynes/cm.) 
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SECTION C 



Finite Length Liquid Column Stability 
in the Presence of a Periodic 
Acceleration Field Oriented Normal 
the Longitudinal Axis 



BACKGROUND -• 


There have been several studies of both forced and free 
oscillations of fluid cylinders. In previous studies which 
incorporated forcing into the problem formulation, the forcing has 
been directed parallel to the longitudinal axis of the fluid 
column. Although those problems are highly relevant, another 
important problem is the case of forcing directed so as to have a 
component perpendicular to the longitudinal axis. 

Therefore, this investigation addresses the stability of the 
fluid column interface in the presence of an acceleration field 
which is directed perpendicular to the longitudinal axis of the 
column. Moreover, this periodic acceleration field gives rise to a 
periodic pressure field as the base state. 


APPROACH 

The analysis is linear and inviscid. Governing equations are 
those of conservation of mass and momentum. The 
nondimensionalization is chosen to be such that the acceleration 
field forcing term will balance the mean pressure field. The mean 
pressure field gradient will then be found to depend on both time 
(t) and theta (the angular dependence) . 

The solution for $ is obtained in the standard fashion. It is 
a solution to Laplace's equation, and since it has been listed in 
previous reports, it will not be repeated here. 

In the solution, there is a dependence of the mean pressure 
field on theta; it is this not possible to consider disturbances 



which are strictly even or odd in theta. Thus, an additional 
summation should be present, with exp(im©) the functional form in 
theta. 

In the case of the finite length column, the velocity normal 
to the end disks must be zero. Also, the anchored triple contact 
line condition is imposed at each of the end disks. 

At the free interface, the kinematic condition and the normal 
force balance must be imposed. The normal force balance will 
contain a term which represents the presence of the forcing, and 
it will multiply the free surface perturbation, which itself is a 
dependent variable. 

Thus, the normal force balance will be non-autonomous due to 
the presence of the sinusoidal forcing function, and a "sin©" 
factor will be introduced, reflecting the mean pressure form. 

The goal is to obtain ordinary differential equations in time. 
The stability of such a system can then be investigated. 

Ideally, the z and © dependencies could be eliminated. For 
example, the z -dependence could be eliminated via introduction of 
a set of expansions in z which would satisfy the anchored triple 
contact line conditions. This is essentially a galerkin method 
approach, and would require use of the orthogonality properties of 
the selected expansion functions. 

The previous difficulties associated with the presence of 
theta have been overcome. Elimination of theta results in an 
infinite system of ordinary differential equations with respect to 
time. Inspection of this infinite system reveals that the even and 



odd azimuthal modes are coupled. 


It was realized during the implementation of these 
considerations that the original basic formulation did not allow 
for the static deformation of the column, independent of any 
wave-like perturbation on the surface. That is, the first 
formulation does not reflect the physical reality as closely as it 
must. 

This was remedied in recent work. The incorporation of the 
static deformation into the formulation is now complete. This 
has had the effect of adding coupling between the longitudinal 
modes as well as the azimuthal modes. 

Floquet theory is then utilized to transform the system of 
differential equations into an infinite set of algebraic equations 
in a parameter, (lower case) lambda. This parameter is the 
eigenvalue. Determination of the eigenvalues will give the 
stability results; if, for all A, Real(X) < 0, then the system is 
stable. The case of marginal stability is given by Real (A) = 0. 
This determination of the eigenvalues must be done numerically. In 
order to effect the numerical investigation, the system must be 
truncated . 

In principle, the finite length and infinite length column 
stability cases could both be done in this manner. However, the 
computer system which is being utilized has restrictions on memory 
and running time which is making the solution of the finite 
problem difficult. 

The aforementioned solution method involves the generation of 



a large matrix representing the truncated system, with an assumed 
Fourier representation in time. The truncations to the solutions 
are done in the summations representing the solution form in 0, z, 
and time, t. For low freguency forcing , the number of terms 
required to for the temporal truncation can be prohibitive. Such a 
matrix can grow so large that it is difficult to utilize an 
eigenvalue solver due to computer limitations. An alternative 
method is sought. 

It is noted that the matrix which is generated is much less 
sparse than it would be if static deformation did not contribute. 

Such an alternative method is in the process of being 
implemented. This involves retaining the derivatives in time; ie, 
assuming no fourier expansion in time, and constructing a 
fundamental matrix which represents the system (itself truncated 
in indices corresponding to 0 and z) . The eigenvalues of this 
fundamental matrix are related to those of the system, allowing 
for calculation of those of the system. Construction of the 
fundamental matrix proceeds via integration of the system for 
suitably chosen sets of independent initial conditions to the end 
of the first forcing period. The solution vector represents a 
column vector in the fundamental matrix. 
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Current interest in microgravity materials processing has focused attention upon the finite 
fluid column. This configuration is used in the modeling of float zones. In this Brief 
Communication, the incompressible inviscid finite length fluid column is subjected to an axial, 
time-dependent disturbance. The response of the fluid system and the resulting interface 
location are determined via Laplace transform methods. 


The presence of vibrations will affect the fluid dynam- 
ical environment in which materials processing occurs. Ex- 
periments performed on Spacelab D1 investigated the ef- 
fect of g-jitter on the stability of finite fluid columns. It was 
found that long columns are particularly sensitive to resid- 
ual accelerations and impulses. 1 The analytical effort has 
been concerned with the axial forcing of slender liquid 
bridges. 2 ' 5 In these analyses, a one-dimensional ‘‘slice’’ 
model has been used. Such a restriction is invalid in the 
situation in which the liquid column is not slender, yet 
whose height and radius are within the static stability lim- 
its for fluid columns. The work of Sanz 6,7 has addressed the 
problem of the unforced modes of oscillation of a nonslen- 
der finite fluid column. 

In this work, the problem of the axially forced finite 
fluid column is investigated, with the axial forcing consid- 
ered periodic in time. The analysis is linear, inviscid, in- 
compressible, and axisymmetric. It utilizes Laplace trans- 
form techniques. The direction of the imposed body force 
is taken parallel to the longitudinal axis, and in the nega- 
tive e, direction. The undisturbed column is cylindrical, 
and is surrounded by an inert gas. 

Quantities will be nondimensionalized as follows: 


panded in terms of Fr about a state of zero mean motion. 
This results in the linearized system of governing equations 

V-u = 0, (3a) 


3(?4_) 

dt 


Vp-g(.t)e r 


(3b) 


with u = V4. For definiteness, g(t) is chosen to be sin(r). 

The boundary conditions at the upper and lower disks 
are that the normal component of the velocity be zero. 
Also, an anchored triple-contact line assumption at both 
disks will further restrict the interface behavior. These are 


34 

dz 



34 

dz 


= 0 , 


A 


(4a) 


/(A,/) =/( — A,r) =0. (4b) 

The constant A acts as a slenderness parameter, and is 
given by (L/2R). As the fluid is incompressible, volume 
must be conserved, and so 



z,i)dz- 0. 


(4c) 


Rx — x, <of 7=f, Ras^—u, paijR^-p, (I) 

and R is the radius of the column. The forcing frequency is 
o/. Pressure and velocity fields are indicated by p and n, 
respectively. The density of the column is given by p. A 
tilde indicates dimensionless quantities. 

Equations are to be rendered nondimensional through 
the use of ( 1 ). This results in 


Because of axisymmetry 



dw 

Tr 


= 0 . 


rmO 


(4d) 


At the fluid column/inert gas interface, both the kine- 
matic and normal force balance conditions must be im- 
posed. These are given by 


V*o = 0, (2a) 

^ + n-Vu = - Vp - Fr gO)?* (2b) 

and the nondimensional interface equation is 

f,=r- l-Fr/(z,r) (2c) 

with f(2,t) representing the functional form of the pertur- 
bation to the interface. The tildes have been dropped for 
convenience. The nondimensional parameter Fr is a 
Froude type number, and is given by (Gq/Ro)}). It is taken 
to be a small parameter. Pressure and velocity are ex- 


^-0. 

(4e) 

Mir) (£4 

(40 


and are evaluated on z=l (A p indicates the change in 
pressure across the interface). The nondimensional param- 
eter B 0 is given by (p/? 3 cj}/y). Gamma denotes the surface 
tension. Note that Eqs. (4e) and (40 have been linearized 
[and are at O(Fr)]. 

Both the system of equations (3) and the boundary/ 
interface conditions (4) are Laplace-transformed with re- 
spect to time. The time dependence is replaced by a pa- 
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rameter s. The initial velocity perturbation is taken to be 
zero. The initial position of the interface is taken to be that 
of the undisturbed cylindrical column /, = 1. This is real- 
istic for the case of a zero mean gravity, which is the case 
in the present investigation. This yields 

jU= -V/»-[l/(^+ l)]e„ (5a) 

<M> 

sF=— . (5b) 

dr 

The upper-case representation indicates the Laplace trans- 
formed variables. 

The velocity potential is found to be 

<30 

<P(W)~ X A„(s)I 0 (k l f)cos[k l ,(z + \)), (6) 

/!= 1 

with k„ = (nir/2\). Here, -4 „(j) indicates the unknown 
coefficients and 7 0 represents the modified Bessel function. 
Use of the kinematic condition yields 


Dim ♦ I “ I U — k\ m „ | ) k} m * |/|(*2m * |) 

— BoJ 0 (k lm . Jr 2 J. 

The conservation of volume requires that 

a(r)[sin( A)/A] + BqK(s) =0. (10a) 

Finally, the anchored triple-contact line condition is 
imposed. This yields a forced set of algebraic equations in 
the coefficients “a(s)" and "6(r),” of the form 

(flll)fl(r) + (al2)Mi)=cl, (10b) 

(a21)a(r) + (o22)b(s)=c2, (10c) 

with all =a21, a22= — al2, and c2= — cl, and with 


sin(A) , (2 sin( A)\ 
(all) =cos(A) — — - — + r 2 f — a — ) 

y Bof 0 (ki m ) 


(lOd) 


00 

sF= X ^,(s)A/,(*,)cos[*,(r + A)J on r=l. (7a) 

m = 1 

However, use of the normal force balance condition yields 
a forced ordinary differential equation for F(zj) at r = l, 

0 +f= ((?Tn) J + W) 

+ X £o*^»U)/o(fr«)coslM*+ A)]. 

»■ 1 

(7b) 


(al2)= -sin(A) + t* ^ 


2 cos(A)^ 
A 

^c/o(^2m+ |) 


) 




(<,,)=, ((?TT)) a | i + (^) 

y ^(/o(*2*+|) 1 

x J : o (Ai. ♦,!>*• + .) ’ 


and with 


(lOe) 


(100 


Note that s appears as a parameter. Also, K(s) is a con- 
stant of integration. The solution of Eq. (7b) is 


F(zj)=a(s)cos(z) + b(s)sin(z) 

* s °*o + (<?TTt) ' 

v (sA n (s)Io{k n )cos[k n (z + A)l\ 

+ M <r^> ) • 

(7c) 


Both forms given for F(zj) must be consistent, and 
must allow for the remaining boundary conditions to be 
satisfied. In order to determine the forms of the coefficients 
“a(s)“ and "b(s)," it is advantageous to expand the func- 
tions sin(z), cos(z), and (z) in terms of cos(/c B (z 
+ A)]. These are found in the literature; 6,7 


A lm = [2 sin(A)/ A]{s/[( l - kljk.j^) 
— B 0 l 0 (k2„)s l ]}a(s), 

1 2 cos(A)' 


( 8 ) 


*2m 


x ( £ ^ 7 r )] (£>! - ,,) "' 


with 


7) f = (* f ( 1 - *J)/,(ik f ) - Vo(* # )^) 

It is found that 


a(s) =0 and b(s) =® (cl/a!2). (lOg) 

Recall that the expression for A lm contained the factor 
a(s). Thus A 2 m will be identically zero, and only the odd 
modes will contribute to the motion. This reflects the na- 
ture of the forcing, i.e., “z” is an odd function with respect 
to the mean plane between the disks. 

The form of the Laplace-transformed perturbation in- 
terface F(zj) is given by 

F(zj) = b(s)sin(z) + B 0 + j - j z 


00 

+ 1 


^2«^i(i)Vo(A2m + i)cos(A:2m + i(« + A)1 


<* ~ k ln,+ \) 


( 11 ) 


with the A lm + | as previously determined. Taking the in- 
verse Laplace transform of Eq. (11) will yield the expres- 
sion for /(z.r). 

In order to gain insight into the form of the perturba- 
tion interface in the presence of a periodic body force di- 
rected along the longitudinal axis of the column, the series 
will be severely truncated. Only one term in each series will 
be retained. The inversion will then be done analytically. 


1629 


Phys. Fluids A, Vol. 3. No. 7. July 1991 


D2 


8rief Communications 1629 




T 0.3 1 .25. 2.25 3.3 


FIG. 1 Deformation fii.t) versus time, for time = 0 25. 1.25, 2.25, 3.25. 
B 0 = 0.50 and A = 0.75. 



Then 


/(z,i)=%in(z)b(t) + 5 0 (z)sin(f) + cos(*,(z + A)] 

x (xTT^T)( 0 X'* (r)si " h| ° ( "' ), '' r 


') 


+ b(0 ] +cos[k } (z + A)) 


m 


X ( J sin(r)sinh[o(f - r) )dr + a sin(r) j , 

V ° 7 (12) 


with 



6(0 -i? -'(*(!)). 

and fP& 1 = a 1 . 


The inverse Laplace transform for F(zj) (see Eq. 
(11)] was taken numerically for different values of (Bo, A). 
Graphical results are presented in Fig. 1. for the case 
(B 0 — 0.5, A = 0.75) at successive Jnondimensional) 
times of 0.2S, 1.25, 2.2$, and 3.25. The maximum value of 
f(z,t ) at each of these times is 0.0019, 0.014S, 0.0131, and 
0.002, respectively. The value of the parameter A indicates 
that the column is not slender. Note the sinusoidal form of 
the interface deformation. This is observed for the case of 
the axially directed forcing. 1 

The equation for F(zj) (Eq. ( 1 1 )] can be viewed as a 
transfer function. Setting s=ia, values of F may be plotted 
in frequency space. This is done, and the results are shown 
in Fig. 2 for A = 2.60, B 0 = 100, and A = 1.25 with 
Bo — 1.00 and 2.00. The curve for the slender column 
(A = 2.60) may be compared with a similar type graph 
found in Ref. 4 (Fig. 3 in that work). It is seen that a spike 
in both graphs occur about the resonant frequency at 
o = 0.34. The long column is shown to be sensitive to 
lower-frequency disturbances (o>< 1), as indicated by the 


FIG. 2. Transfer function F(z*.u). *• = - A/2, for the cases (A = 2 •>'. 
8,= I 00). (A = 1.25. 8 0 = 1.00), and (A = 1.25. B 0 = 2.00) 


relatively large functional values (for F(z,j)] in this work 
and in Ref. 4. Note that the numerical value of F(z*ju) 
(with z* — — A/2) compares well with that of the defined 
transfer function perturbation amplitude in Ref. 4. 

It was found that nonslender columns are not similar!) 
sensitive to these low-frequency disturbances. This is 
shown graphically in Fig. 2 for A = 1.25, B 0 = 1.00, and 
B 0 — 2.00. Note the low numerical values of F(z*,u) as 
compared with that of the case A = 2.60, B 0 = 1.00. 

The work of Sanz 4 determined the nondimensiona^ 
natural frequency a u , of the fluid column for a range of A 
values. In that work, time was nondimensionalized b> 
(plO/y) 1/5 - It can be shown that this results in the follow 
ing equivalence (Of/utl^) = B 0 ( 1/5^,). 

The perturbation interface was determined (numeri 
cally) in a number of cases for which the (Bo,A) value 
corresponded to eigenvalues of the unforced system. Ir 
these cases, the inviscid finite length fluid column wa 
found to be very sensitive to forcing. 

Various orders of magnitude for the vibration leve 
have appeared in the literature.* The range of interest o' 
parameter B 0 can be constructed for fluid cylinders of i 
given density and radius. The interface shape can be ob 
tained for specific values of B 0 and A. 
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Fig. 2. Typical CD cocffldeat tiaw histories (AF * 7J pd/91.7 kfs; 
dcu; without k«), 



Fl|. 3, Depcodetct of drag coefficient oa dmtloi UkP* 7.5 pri/ 
91.7 fcJTm; deaa; without berm). 



elostion where base scavenging due lo vortex shedding is 
substantially reduced {HE/H « 0.25. Fig. 2a). 

Because of the sensitivity of base pressure to the vortex 
shedding phenomena, investigators have pointed out that re- 
duction in base drag can be achieved by modest geometric 
changes that interfere with the vortex generating mechanism. 
Such configuration modifications include fineness ratio, 1 cor* 
ner rounding,* splitter-plate length, 5-5 trailing-edge spoiler, 5 
and surface proximity (current results). Note, for example, the 
onset at Late time of periodicity in the CDB signal at HE/ 
H = 6 in Fig. 2b (the vortex shedding “signature*’) and the 
absence of any periodicity for the data for HE/H = 0.25 (e g.. 
Fig. 2a). Since the drag “beat** frequency is typically twice 
that of the comer shedding frequency, 10 an estimate can be 
made as to the Strouhal number (5) for the current study 
S = *H/u m 

where h represents the Strouhal frequency (»0.5x drag peri- 
odicity), H the model height, and u m the local freestream 
velocity. The computed value so derived is seen to be some- 
what higher (0.15) than previous measurements (0.13) (Ref. 
3). The onset of CDB periodicity is also consistent with the 
observed increase in base drag as would be expected since base 
scavenging is controlled by the establishment of the von 
KirmAn vortex street. 5 4 
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Fluid Column Stability in the Presence 
of Periodic Accelerations 


critical to loads evaluation studies as fineness ratio. Reynolds 
number, and comer rounding considerations. Note the favor* 
able comparison between the current steady-state data for the 
two elevation limit cases, HE-0 and HE-*, with earlier 
results.*'*' 1 * Apparently the influence of boundary-layer 
buildup ( 6 , 0 *, * 0.3 in./0.0l27 m) on measured CD results is 
“small” as demonstrated by the favorable agreement of the 
low HE data with established splitter-plate results. Also, the 
effects of vortex shedding (evidenced by the presence or ab- 
sence of periodicity in the CDB data), as already discussed and 
reviewed by Bearman* and Vickery," can be seen to play a 
major role in maintaining low base pressures (high CDB) at 
high elevation ( HE/H = 6 , Fig. 2b) and low CDB at low 
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Introduction 

T HE float zone configuration is used in crystal growth. It 
may be modeled as a liquid column held by surface ten- 
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sion forces bctwe<n two end disks. In the case of crystal 
growth, the thermal and soIutaJ fields as well as those of 
velocity and pressure arc needed to characteri re the physics. A 
compelling reason for crystal growth experiments in a micro- 
gravity environment such as that onboard the Space Shuttle is 
that buoyancy forces are greatly reduced. 1 However, this envi- 
ronment is not quiescent due to the presence of impulse type 
disturbances from small thruster firings as well as those from 
periodic vibrations. Given certain * , environmental ,, condi- 
tions, such as the presence of a periodic acceleration field, it is 
possible that a fluid dynamical instability would develop. This 
would adversely effect the crystal growth. 

The crystal growth environment cannot be separated from 
the fluid dynamics of the liquid column. This latter topic is 
the focus of the present work. The question of column inter- 
face stability in the presence of a periodic acceleration field 
having a component normal to the longitudinal axis of the 
isothermal cylinder is investigated. The fluid column is taken 
to be infinite in length. Roquet theory is used in the stability 
investigation. 

Previous work has determined the natural oscillations of the 
liquid column. This work was done first for the case of the 
infinite column 2 and later for the finite length case. In the 
latter work, both axisymmetric and nonaxisymmetric oscilla- 
tions were considered. 34 Results for the infinite length case 
were found to be good approximations to those for the finite 
length column, both numerically and with regard to trends. 

Interface behavior of the finite length liquid column in the 
presence of time-dependent forcing has been investigated for 
tlie case in which the forcing was parallel to the longitudinal 
axis of the column. 34 The investigations considered the 
column behavior subject to a sin(/) forcing for both the invis- 
cid case of general aspect ratio (within static stability limits 3 ) 
and the viscous case in the slender column limit. 4 

The problem of interface stability of the fluid column in the 
presence of a periodic acceleration field that has a component 
normal to the longitudinal axis of the column has not been 
investigated. Previous work has considered the interface sta- 
bility of a highly idealized infinite slab-like configuration in 
the presence of a periodic acceleration field oriented normal to 
the interface. 7 Interestingly, this study was motivated by ex- 
periments in microgravity. 

Use of the infinite length configuration in this stability 
stndy results in a simplification in that the standard boundary 
conditions at the solid end disks are not applied, and the focus 
remains on the interface stability. If an extension of this work 
to the finite length configuration is of interest. Roquet analy- 
sis would be appropriate, although the implementation would 
be more complicated. 


Formulation 

The basic configuration is that of an infinite fluid column of 
circular cross section. The fluid is incompressible, and the 
surrounding medium is of negligible density. Perturbations 
are taken to be irrotational. The analysis is linear and invisdd, 
with the nondimensionalized governing equations those of 
continuity and conservation of momentum (linearized Euler). 

The frequency of the periodic forcing is denoted by «/. 
Pressure and velocity fields are given by p and «i, nondimen- 
siooalized as follows: 

Ri ~ x w/*r» t Rufi~u p(Ru/?fl=p (1) 


The mean state is one of zero velocity; however, the mean 
pressure is time dependent. Consider a wave-like perturbation 
propagating on the interface. The governing equations for the 
perturbation are developed as follows. Expand in the small 
perturbation parameter i 


+ (3) 

Substitution of Eq. (3) into Eqs. (2) yields the mean system 


= -Fr cosl/)Vl(l - r)s.n$l (4) 

(with the parameter Fr of order one) and the (order e) pertur- 
bation equations 


V 2 * = 0 


(5a) 


3/ 


Vp 


(5b) 


The spatial dependence of the velocity potential can be deter- 
mined via solution of Eq. (5a), which is 


Hr, 3, z, t) = E^(0/^(^)cxp(^)cxp0>n3) 


( 6 ) 


f m is the m th modified Bessel function. Clearly, the solution 
involves a superposition of the azimuthal modes. 

It is through the boundary conditions that the free interface 
can be determined and the stability characteristics investi- 
gated. Let the equilibrium interface be given by 

Ft » r - 1 - «t(0, z, t) • r - 1 - <EC(/)exp{/*z 4- im$) (7) 


The kinematic condition and the normal force balance at the 
free interface must be satisfied. In addition, the requirement 
of conservation of mass, which reduces to a conservation of 
volume condition, must hold. 

The linearized kinematic condition (at order <) is 

(8a) 

ot 

at r = 1 with u, » d<t>/dr. This results in 

^^exp(/*z + imt) • El/„(*)l'«4(f)exp(/** + imt) (8b) 

The normal force balance requires that the difference in 
pressure across the interface be equal to the curvature multi- 
plied by the surface tension force. In nondimensional form, 
this is given by 

Bo x A( p- t - n , + tp) » V ■ ■ (9a) 

with ■ * VFe/tvFel. Bo is the nondimensional parameter 
with i the surface tension and p the density. At 
order «, the linearized form of this interfacial condition can be 
expressed as 

1 + lo + f# + IBo cos (/) sin djij * (9W 

Fr, required to be of order one. has been set equal to unity- 
The subscripts indicate partial differentiation. This yields 

El(l - m 1 - k*) + Bo co $</) sin *| x C(/)exp(/*z + 



Tildes indicate nondimensional quantities. 

The continuity and Euler equations are then 

V-fi = 0 (2a) 

vfi- - vp -f>cos(/)v((l -flsind) (2b) 

ot 

with Fr - (Go/RJ,) a Froude type number. Go is the ampli- 
tude of the periodic acceleration field. The functional form of 
the time-dependent forcing is selected to be cos(f). 

06 


- EA (f)/-(*)exp(i*z + imt) 

The sin 9 dependence can be re-expressed as an exponential 
function. Then 

C„(/X1 - m* - **) + (y)(-0 cos(r)[C._ ,(/) 

-C..,(/)| = *oU*)^=) (9d) 
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a nd azimuthal mode coupling occurs. Use of Eqs. (8b) and 
(9d) yields a nonautonomous second-order equation in C„(f), 
«ith mode coupling (^.indicated by the subscripts m ) 

Cl - 10 -m* -k 2 )/ 8o\\ll/ ! m \C m 

- I H/2t m )cos (fK - i)|C. . ,(/) - c m . ,(/)) (10) 

Keep in mind that cos(f) can be rewritten in exponential form. 

]t is at this juncture that Roquet analysis is applied. For 
convenience, let £„(/) = [dC„(f)'df)- Then take 

[C„(f). £„(/)) = E;:MC.,, Em.i) * expl(X ♦ il)t\ (II) 

The constant coefficients C mJ and E mJ are unknown. The 
Floquet exponent is denoted by the eigenvalue X, which is in 
general a complex number, and which is also unknown. The 
nonautonomous differential system is thus transformed into a 
homogeneous algebraic system for (C mJ , E ml ) with the un- 
known parameter (eigenvalue) X. If Real(X) is greater than 
zero, the interface of the cylindrical column is unstable to the 
growing wave-like perturbation. Of course, the milieu in 
which this disturbance is propagating includes the periodic 
base state pressure. 

Use of Eq. (11) in Eq. (10) (rewritten as two first-order 
modes) yields the infinite algebraic system given by 

(X + il)C m j = E m .i (12a) 

(X + //)£«., = (d - «* - 4rty(A»l</;/ 

+ o/«/;/i.x - /xc._ , ♦ c m . u . , 

-C,*u-i"Ci*u*i) (12b) 

Note that the harmonic modes (indicated by I) as well as the 
azimuthal inodes (indicated by m) are coupled to both their 
preceding and successive modes. 

Several remarks are in order concerning the truncation. 
Once the truncation in m is done, the number of azimuthal 
modes that contribute are fixed. It is to this truncated system 
that Floquet analysis is actually applied. To obtain numerical 
values for X, it is necessary to truncate the number of harmonic 
modes in time, i.e., the range of / values. The eigenvalue 
problem is, therefore, a problem of the truncated system. 

Results 

The results pertain to the eigenvalue solutions of system 
(12a) and (12b). NAG library routines were used in determin- 
ing the eigenvalues. Truncation values of L - 1151, that is, 
- IS s / s 13, and M » 14 were found to be sufficient. Wave 
number values ranged from k «0.10 to * » 3.00. The parame- 
ter Bo was varied from 0.01 to 10.00. 

Forftct .0, the interface is unstable to the wave-like pertur- 
bation in the presence of e mean periodic acceleration field 
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over the range of Bo values considered. (Note that this range 
encompasses four orders of magnitude; see Fig. 1.) 

As k is increased, such that 1.00 s 4r s 1.20, the interface 
remains unstable to the perturbation for the two larger Bo 
values, equal to 10.00 and 1.00. However, marginal stability 
(Real(X) » OJ ensues at the smaller Bo values. A decrease in Bo 
can be interpreted physically as an increase in the surface 
tension. So. as the surface tension increases, the restoring 
force is sufficient to result in marginal stability over this range 
of wave numbers. 

Perturbations corresponding to the larger wave numbers 
(smaller wavelengths) that were considered. 2.0 sk s 3.0, 
were found not to grow in time for Bo = 0. 01 -1.00. That is, 
instability [with Real(X)>0) of the interface to perturbations 
of these larger wave numbers (and smaller wavelengths) occurs 
only for Bo = 10.00; otherwise, Real(X) = 0. An alternative 
physical interpretation to that involving a variation in surface 
tension for differing Bo values can be developed. Since Fr was 
taken to be unity, wj (the forcing frequency) is proportional to 
Co, the amplitude of the periodic acceleration field. Utilizing 
this relation in the definition for Bo yields Boa(Go/y). For 
fixed surface tension (and, of course, density and column 
radius) values, an increase in Bo would result from an increase 
in forcing amplitude. It is at the highest such amplitude con- 
sidered that the interface was found to be unstable [with 
Rea!(X)>0] to (he perturbation. 

It is noted that the range of Bo values used corresponds to 
values of Go and wy that would be of interest in a microgravity 
environment for certain ranges of surface tension values. 
(Roughly, Go s 10 and 0.5 Hz<«y<5Hz 

for y values of 1-100 dynes/cm.) 
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Coherent Structure Interactions 
in Excited Coaxial Jet of 
Mean Velocity Ratio of 0.3 
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Nomenclature 

D = outer diameter of outer nozzle 

d * inner diameter of inner nozzle 


Fig. 1 Stability diagram: the stability of the configuration to wavo- 
Hlte pertorhatious for a range of Bo parameter values n wave lumber 
M b shown. The cross-hatched area Indicates uo stable regions la which 
disturbances arc growing la time. The remaining area corresponds to 
that of (he marginal stability state. 


Received April 7, 1992; revision received Dec. 21, 1992; accepted for 
publication Dec. 28, 1992. Copyright €> 1993 by the American 
Institute of Aeronautics and Astronautics, Inc. All rights reserved. 
'Graduate Student, Department of Mechanical Engineering. 
fProfessor, Department of Mechanical Engineering. 


07 


; 


AlAA JOURNAL 


* OL 19. NO i i 


Effect of Periodic Accelerations on Interface Stability in a 
Multilayered Fluid Configuration 


M. J. Lyell* and Michael Roht 
West Virginia University, Morgantown, West Virginia 26505 


The locreislng number of research opportunities In a mlcrograrily environment will benefit not only fun- 
damental studies tn fluid dynamics, but also technological applications such as those Involving materials pro* 
cessing. U particular, fluid configurations that Involve fluid-fluid Interfaces moetd occur In a variety of exper- 
imental Investigations. This work Investigates the stability of a configuration Involving fluid-fluid interfaces tn 
the presence of a Ume -dependent (periodic) forcing. The fluid configuration Is mutitbv ered and infinite in extent. 
The ana}) sis b linear and Invlstld, and the acceleration vector b oriented perpendicular to each interface. A 
Roquet analysts b employed, and the resulting algebraic elgensystem b truncated. Nondime rulonal parameters 
appear b the algebraic system. A numerical study b performed to elucidate the regions of Instability and the 
effect of parameter variation on the fluid configuration stability. 


Nomenclature 

= matrices representing final system 
A(f) = time -dependent coefficient, nondimensional 
Bo 8 nondimensional parameter, ** fa^JFGJy) 

Fe = equilibrium interface, nondimensionalized 

Fr 8 nondimensional parameter, « \GJ{Huj)] 

C 0 8 peak value of forcing, dimensional 

g(/) 8 forcing term, nondimensionalized 

H 8 height of middle slab, dimensional 

k 8 wave number of perturbation, nondimensional 

A = outward pointing unit normal to interface 

p = pressure field, nondimensionalized 

ir = velocity field, nondimensionalized 

I = time, nondimensionalized 

y = surface tension 

X 8 Floquet exponent, eigenvalue 

p = density 

4 = velodty potential, nondimensionalized 

<jy = frequency of periodic acceleration (forcing) 

Subscripts 

1 8 (finite) middle layer fluid region of height H 

II 8 upper layer (unbounded) fluid region 

III 8 lower layer (unbounded) fluid region 

2 8 upper interface 

3 8 lower interface 

Superscripts 

() 8 differentiation 

() = unit vector (e.g., A or part of vector defi- 

nition (e.g.,1) 

O 8 dimensional quantities 
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Introduction 

C URRENT interest in microgravity materials processing 
has focused attention upon certain relevant aspects of 
fluid mechanics in this environment. In particular, a number 
of materials processing applications involve fluid-fluid inter- 
faces. 

The environment onboard the Space Shuttle is not strictly 
a microgravity environment. Rather, residual accelerations 
exist that could affect any ongoing materials science or space 
processing experiments. A recent summary 1 indicates that 
accelerations indude those in the frequency range of <10 Hz, 
with acceleration levels ranging from 10’*g o# * to 
In addition to periodic forcing, residual accelerations may be 
of impulse type, due to such causes as stationkeeping ma- 
neuvers and astronaut motion. 1 

This work investigates the effect of periodic accelerations 
on the interface stablity of an idealized fluid configuration. 
The fluid configuration is multilayered and infinite in extent 
(see Ftg. 1). The accelerations are periodic about a zero mean 
g level and are oriented normal (in the i t direction) to each 
interface. 

Previous wort has investigated the stability of a single pla nar 
free surface subject to periodic forcing in the direction per- 
pendicular to the interface. 14 These studies were both done 
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u.ncf. In t*c c? 3:-..s-ii i-j (.'tv?:!.* *e c^ i -er 
«« cylmd'icj" *‘:>c T>f.r *- 3 :.*:% ?c-j‘»ed.-. a Mi* 1 ?u 
equation that £;%c:-fi Jfse ?.me deponent a-;* :*ie cl t*€ 
disturbance t' ;> *?re aMe »o r^aVc sutcffev* «w *j 
the interface :y ti*ed on Ijv*« rtather areal ;tcp- 
c mes of Mjth.ea flattens. The case of a reetar.gvi'ar cots* 
uiner has been addressed recently.* and the results extended 
into the nonlinear regime. Both of these investigations united 
an inviscid analysis 

The effect of viscosity on the stability properties has been 
studied recently in tdealned infinite Of semi-infinite config- 
orations that had one fluid-fluid interface.’* Again. the pe- 
riodic forcing was directed normal to the interface. One the- 
oretical investigation was done assuming a rero mean g level 
and pertains to the microgravily environment.* Results mere 
obtained for different repons of parameter space, and tome 
stability boundaries were determined. 

The present work wifl study the stability of an idealitcd 
multiple-layer fluid configuration that is infinite in the i $ and 
i directions (as well as unbounded in the i a direction) and 
is subjected to periodic forcing ooroial to each interface. The 
height of the middle layer is finite, to general, the densities 
in each region differ, as do the surface tension values at the 
upper and lower interfaces. The analysis is linear and invisrid. 
A Roquet analysis is employed. Details of the analysis are 
presented in the nett section. 

The configuration is idealized. as was the case in previous 
work. However, the use of oondime nsional parameters ia the 
stability investigation wffl enable treods to be Accrued and 
parameter regions (of instability) to be identified. 

Equation Development 

Governing Equathm 

The to v« mint equations are those of conservation of mass 
and momentum applied to an invnbd and incompressible 
fluid, nod are 

f • « - 0 0) 

♦ I'^lj ■ (1) 

The time -dependent body force tens is indicated in Eq. (2). 
C. represents the peak value of the acceleration due to the 
periodic forcsnf. The hmetioo {(•/) •» periotfe a«l nil be 
takea to be e cosine function. The fluid system • fineanted 
about e state of zero mean motion, and quadnticalty small 
terms (after espansioo ia a small parameter c) are neglected. 
Use of nondlmeosionafizJtiom 


fHw. 

i - ny*’ 1 

(3*) 

t » Hm/t, 

fpJiWf 

(») 

V • a - 0 

(4) 


;)a:e\;*ed i ^e - li e -r;*r a-* *: .cr reg . to >»e?d 

<6. = {4.:lr%r(t:) * 5.;)fip( -l:::evp-^ • -,J) (6s) 

< 6.1 ° iCfr) r%p( •4;*| evp* ♦ •*»*?) ( 6 b) 

d„, * (D(r)eip(l.*iJevp(fiif ♦ mr]) (6c) 

Note that A. 0. C. and 0 are time dependent coefficients. 
The w*\e number k is gi*en by v(£ ♦«?). 

Equations representing the upper and tower equilibrium 
interfaces arc 

fir, » I - 1 - ( £(#) cip(ij(r ♦ my)) (7) 

/>,«*-« f\i) cspOJlr ♦ ^yj) (g) 

Note that £ and Tare time -dependent coefficients. Thus, for 
example, the functional form of the perturbation to the oeaa 
interface at r - 0 is given by the expression £(i) e xp(j{6r ♦ 
my)). At this stage, the time-dependent coefficients are un- 
knowns. 

Boundary Condition* 

The kinematic condition, which holds at each interface, 
requires that 

♦ ■ • Vfc - 0 (9a) 

After linearization, only the ». velocity component win coo- 
tribute. This results ia 

£(t) ♦ k C{i) e-* - 0 oo r - I (9b) 

-t{t) ♦ *0(0-0 oo « -• (9c) 

In addition to the kinematic condtfon, the normal com- 
ponent of the velocity is continuous across an interface, which 
yields 

(*) - ( **) * * • 1 < 10 *> 

($) ■ (i?) * (,0b) 

or 

A(i)t* - 0(f)- -C(i) oo Ml ; 

AM - 0(r) - D(i) on a • • i. 

The remaining boundary condition is the (finearhed) normal 
force balance across an interface. In aondimemaona! form, it 
is 

(0o/fr)(dp) - V-4 (II) 

The unit vector 4 ts the frnearued outward poiting normal 
to the interface. After substitution, Eq. (11) yidft 


Tbe parameter Ft - (PJH*, 1 ) a tabeajobe ro^My of order 
1 Note that the mein pressure field win be periodic in feme. 
fAlso , is expanded into both a mcaa and pettvbatioa coo- 
JS£ri£) A potential kinetic ♦ with . - ?*. is defimd. 

Substitution into Eq. W y*** •&****» 

bations arc taken to be wavelOe in the (zy) plane. That a. 
they are osoltalOty. The tesultin* dMlerential equatwm (in Q 
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♦ - (FriBoJVEW on r-1 

(12n) 

(Op. - PiuKpJ/>*(OfW ♦ <P/Pjtftt ♦ *£ 

-(pu/P-)^(i))-(/WfloJk*fW c z-o 
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Bo 2 = (p.. If Cj^i) 


Bo , = (p.. tf GJ 1} ) 


Fr = GJ(H», : ) 

The term y.(yj) indicates the surface tension at the upper 
(lower) interface. 


Final Sjstem of Equations 

Through the utilization of Eqs. (10c) and (lOd), the C(f) 
and D(t) coefficients can be eliminated to yield the following 
system: 

A(/){-3(l + Pj,)/( 1 + P„ + P,.)}f* 

+ B(/){3(p., - J)/(l + Pj. ♦ P,.))*'* 

+ £(0 |3(P;i - l)/(l + P„ + P,|)} Frg(l) 

= ( FrlBo : ) k'- £(/) 

A{i ) {3(1 - p„V(l 4 p„ + P,,)l 

♦ ^(0 {3(1 + P»iV(l + Pai + P>i)l 

4 F[ i) {3(1 - p„)/(l + p„ + Pm)) Frg(l) 

= (FrlBoJ k> F[l) 

Fit) - * MW «• - *(')<**} 

Fit) - k[A(t) - £(/)} 

with p,, = (p,/p,) and p„ « (Pu/Pi)- 

The governing system of equations has been reduced to 
four ordinary differential equations (in time) with noncoo* 
slant coefficients and with four unknowns. It remains to de- 
termine the stability of this system for various values of the 
parameters Fr, Be„ Bo y , (hi. p>i. and k. As was mentioned 
previously, g(r) is chosen to be cos(r). 

Floquet theory’ can be applied to system (13). Let the time- 
dependent coefficients be expressed as 

{Ail).Bil),Eii).m - £ {A„B..E„Fj €~€" (14) 

• • -• 

The system of four ordinary differential equations in time 
is now re-expressed as an infinite algebraic system, with the 
Floquet exponent appearing as another parameter. That is. 


(13a) 


(13b) 

(13c) 

(13d) 


<*(X + «){- 3(1 + fc,yu + (hi ♦ fc,)M. 


+ *‘(X + <n){3(pi, - iy(l + p,, + P„))fl. 
- (Fr/BoJk'E. 4 {3(pj, - iy(l 4 p,, 

(15a) 

(X 4 «){3(1-P„y(i + Pll ♦?„)}*. 


4 (X 4 in) {in 4 p„y(l 4- p,, 4 p,,)} B. 

- iFrIBoJk’F. 4 {3(1 - p„y(l 4 Pjl 4 p„)) 
x (Fr/Tj{F m - 1 4 F 1 , , ,) ** 0 

(15b) 

(X 4 in) E. 4 k e* 1 B m - k t* A. » 0 

(15c) 

(X 4 in)F m 4 kB.- kA, = 0 

0*0 
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The value cf n varies from -* to The infinite set o' 
Komcgcncc-s equations can be vxrittcn in matrix form as 

l£ l (4) - x [\i = 0 (Ua> 

or 

= ( 16 k) 

The lalter representation is in the form of a generalized ei- 
genvalue. problem. The Floquet exponent X serves as the ei- 
genvalue. The infinite column vector t contains (A m ,B mt E mt FJ T . 
For details concerning the form of matrices A and £, and 
column vector i, see the Appendix. 

The sign on the Floquet cxponcnl/eigenvalue X indicates 
the stability. If X is a positive value, then the disturbance will 
grow and the interface configuration will be unstable. As this 
is a linear analysis, no information can be obtained concerning 
the finite amplitude (nonlinear) form of the configuration. 
Numerical methods are used to determine the value of X for 
different values of the parameters. These arc discussed in the 
next section along with various results. 

Results 

t Comments oo Numerical Method 

The linear algebraic system given by Eqs. (16) is truncated. 
Results are presented for a truncation value of n « |25|, whkfc 
involves a system of 204 equations. This resulting system forms 
a generalized eigenvalue problem, with the matrices A *nd 
fl having complex entries. At this level of truncation, both 
A and £ are sparse. 

Eigenvalue techniques were used to determine X values. 
An IMSL routine (DGVLCG), based on an LZ algorithm, 1 
was utilized. In particular, we are interested in the value of 
the largest ReaJ(X), which gives the fastest growing mode. If 
Real(X) is positive, the perturbation is growing exponentially. 
Hence, the fluid configuration is then unstable with respect 
to the perturbations in the presence of the unsteady periodfc 
acceleration (forcing) field. 

Several checks were made to insure that the eigenvalues 
obtained were correct. Higher truncations yielded the same 
largest Real(X) value. As an additional check, the generalized 
eigenvalue problem was reformulated as a standard eigen- 
value problem of the form (£~ ( 4 - X /) J = • and the ei- 
genvalues of this new system were determined. This was done 
by utilizing an alternate IMSL routine (DEVLCG) based cm 
a different numerical algorithm. Again, the values of ReaJ(X) 
were in agreement with those obtained previously. As a final 
check, a limit .case was run (for various parameter values). 
The value of Bo t was set to infinity, which is equivalent lo 
setting the surface tension at the upper interface to zero. In 
addition, the densities in regions I and II were set to the same 
value. This limit case represents the one interface case. The 
actual one interface case was developed separately. Results 
obtained from using the two interface system (and code de- 
veloped foe it) in the aforementioned limit case and from 
solving directly the one interface system agreed quite closely. 

Dtscusstoo of Results 

A parametric study was performed to investigate (he effects 
of variation in Bo a , Bo> % Fr % and the density ratios with wave 
number. Because of space limitations, a subset of results that 
were obtained are presented. (These are illustrative and typ- 
ical.) In each case, if the periodic forcing function g(r) were 
to be set to zero, and with the mean gravity level zero, the 
configuration would be stable to the wave-type small ampli- 
tude disturbances. That is, the interface would simply oscil- 
late. It b only with the forcing, and in the indicated parameter 
regions, that instability does result. 

The values of the wave numbers vary between 0.10 and 
5.00. Although not all results are presented, the investigation 
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considered ffoj values of (1 .0, 0.1, 0.01, 0.001). The values 
of Bo t were set equal to Bo it or were twice those of Bo l in 
different studies. 

Figures 2-10 show graphs of Real(X) vs wave number for 
various values of Soj, Bo t , Fr, fa, and p,,. Note that Real(X) 
is the largest eigenvalue and represents the fastest growing 
(unstable) wave if positive. Regions for which Real(X) £ 0 
are indicated by being set equal to tero in the graphs. That 
is, the value of Real(X) < 0 is not of interest. 

In Figs. 2-7, the effect of Bo 3 and Bo 3 oo stability for 
different Fr and density ratos is elucidated. In Fig. 2. Bo, is 
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Fig. 7 E/Tect of Bo oo stabtrity: p u « 1.00; p M « 0.50; /> = 1.00; 
Joj a 2 x (JoJ. 

set equal to Bo t for tbe density ratios fa » 1.0 and fa ■ 
0.001225 and for Fr = 5.00. the unstable wave number region 
is broadest for the largest Bo^Bo^ values. As Bo^BoJ is 
decreased, the unstable wave number region shrinks to en- 
compass fewer k values and tends toward the lower k region. 
For low Bo^BoJ values, this two interface configuration is 
unstable to the longer wavelength disturbances in the presence 
of periodic forcing. The effect of a decrease in Fr, keeping 
the other parameter values the same, is seen in Fig. 3. Tbe 
range of unstable wave numberi broadens, with smaller wave- 
lengths falling into the unstable region. 
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The density ratios indicated in Figs 4 and 5 pertain to a 
gas !iqutd.gas configuration As was the situation in the pre- 
ceding two figures, the broadest range of unstable wa\c num- 
ber values occurs for the largest 80 .( 80 ,) value. As Ft is 
decreased, the unstable w3ve number band encompasses larger 
values (corresponding to smaller wavelength disturbances). 
Inspection of Fig. 5 shows a slight separation in unstable wave 
number bands near k = 2.175. This k value falls into an 
unstable region at the smaller Ft value given in Fig. 4. Note 
that the configuration represented in Figs. 4 and 5 is stable 
when Bo : (Bo } ) is 0.001. That is, Real(X) > 0 only for the 
larger 805 ( 80 ,) values in the indicated k ranges. 

The effect of Bo 2 not equal to Bo, on the configuration 
stability is shown in Figs. 6 and 7 (for the given values of Fr 
and density ratios). In F»g. 6 . the Bo values are equal. This 
is changed in Fig. 7, with Bo* twice Bo ; . Physically, the in- 
crease of Bo, w hile keeping Fr fixed can be interpreted as a 
change (decrease) in the surface tension value at the lower 
interface. The predominant effect is to broaden the range of 
unstable wave numbers for each set of Bo : , Bo, values. Note 
that the numerical value of the real part of the Floquet ex- 
ponent X is increased for Bo, twice the value of Bo ; . indicating 
a faster growing “fastest growing’* disturbance at the lower 
Bo : values. 

In Fig. 8 , the effect of holding the (Bo ; .Bo,) values fixed 
(for the specified density ratios) and varying Fr is shown. As 
Fr is decreased, the range of unstable wave numbers increases. 
Physically, this can be interpreted as a decrease in configu- 
ration stability with respect to the wavelike disturbances for 
somewhat larger frequencies of the periodic forcing. 



Fig. $ Effect of Fr on stability: p u * f.00; p u = 0.001225; Bo, « 
0.01; Bo, s 0.005* 
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Fig. 8 Effect of height variation: Fr « 5.8 cm-H; Bo, » Bo, » 
0.001% x H’cm-*; p* * 1.80; p M * 1.28. 


The height of the finite middle slab a i physical quarvy 
that appears in both the Bo and Fr r.c* dime ns ion a I param- 
eters. In particular, Fr is inversely proportional to the hc; w .t. 
whereas Bo depends on the square of the height. An increi^e 
in height H implies a decrease in Fr and an increase in Bo. 
As was seen in Figs. 2-7, the region cf instability in wave 
number space becomes larger as Bo increases. From Fig 8 . 
it is seen that the broadest region of instability (with respect 
to k) corresponds to the smallest Fr values. It is anticipated 
that an increase in H will result in an unstable configuration 
for a broader range of k values. This is borne out in Fig. 9. 
Results are presented graphically for the case G # = (10 " 4 
£e«s*)« «/ 3 0.1 Hi. andy, « y, = 50dynes/cm. In addition, 
fct « 0.8 and p,! * 1 . 2 . 

The effect of density ratio difference on stability is pre- 
sented in Fig. 10. Values of p J( and p„ represent the density 
ratios of the upper and lower regions to (hat of the middle 
layer, respectively. Among cases indicated, the largest mag- 
nitude difference in the density ratios is (p, 4 - P:* = 9). This 
corresponds to the case having the largest range of unstable 
wave numbers. Note that all three cases belong to a family 
with p,, » 1.0. Id addition, the case io which both density 
ratios were set to I, indicating equal densities in all three 
regions, was addressed. Under the action of periodic forcing, 
lack of density differences among the layers resalts in lack of 
instability. 

The wave number at which the subbannonic occurs is plot- 
ted in Fig. 11 for a range of Fr values at the indicated Bo 
values and density ratios. It is seen that there b a shift of the 
subbannonic to lower wave numbers as Fr increases, i.e., as 
the forcing frequency decreases. 



*•100 *-100 * 1 * 1.00 

*•10.00 *4.50 *4.001225 

Fig. 10 Effect of deoslly ratios on stablfit): Bo, a M0 x p„; Bo, 
» 0.05 X p„; Fr * 1.00. 


Subharmonk Case 



Fig. It SubhumooJc mode at (Fr,*) valves for Bo, a Bo, « 0.01: 
Pu ® 1-00; p M a 0.001125. 
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F*. 12 Stability boundaries: cases p u * p M 3 0.001225; and ^ 
*= 1.00, Pil = 0.001225. 
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Fig. 13 Subifity boundaries: cases fa = 1.00, fa « 10.00; aod fa 
° 100, p M = 0.001225* 


low A region. This conclus on '‘as determined through con- 
sideration of a limit ease that imol'cd setting Bo> to » r. T: r. . : y 
and p,, to 1 in order to studs the effect of having onl) one 
interface using the code developed for the two interface con- 
figuration. In this way, the parameters are consistent in both 
problems. 

In the stability boundary diagrams of Figs. 12 and 13. the 
rectangular filled regions indicate instability. It is clear that, 
over the indicated range of Bo : (BoJ and Fr values, the con- 
figuration is generally unstable at lower wave numbers {k < 
1). An exception is the gaj-liquid-gas configuration (fa = fa 
* 0.001225) at lower Fr values. Note also that the regions of 
instability are punctuated by stable regions for the range of 
Bo values and ail but the highest Fr value. 

The configuration used in this study was idealized. In m 
actual space processing application, the fluid configuration 
would not be infinite io extent. Boundary conditions pertinent 
to the specific application mould have to be taken into ac- 
count. Nevertheless, the results obtained in this study can be 
used in a qualitative manner when considering a specific ma- 
terials processing geometry. 

For example, it has been determined that the multilayered 
fluid system is, in genera), unstable over a broader range of 
k values than the one interface fluid configuration. This has 
implications for a float xooe processing technique in which 
the fluid cylinder is multilayered. Abo, the investigation into 
the important subharmonic case shows that it occurs at wave 
numbers that increase with decreasing Fr values. That is, io 
the more unstable Fr range, the subharmonic (Floquet ex* 
pooent X = 1) will occur at smaller values of the perturbation 
wavelength. A materiab processing fluid configuration then 
could be susceptible to instability due to small wavelength 
fluctuations in the presence of periodic forcing. The nondi* 
mensionai idealized system has been studied over the range 
of parameter values relevant to a microgravity environment 
(including the Space Shuttle), as can be seen bom the Intro- 
duction. Configurations involving fluids of specific interest 
can be investigated at greater length using this methodology. 

Appendix 

The form of the column vector 1 is given by 


Stability boundaries of Fr vs k are plotted in Figs. 12 and 
13. (Recall that Fr is inversely proportional to the square of 
*y) This is done for configurations of different density ratios 
(as indicated by the different area fill patterns). Moreover, 
on each graph, multiple values of Bo Zl Bo % are represented. 
The unstable regions are indicated by the rectangular filled 
regions. No meaning is ascribed to the width of the rectangles. 


* - (- - • Y CA1) 

The form of matrix is given by 


-..-1 0 0 0 

...0-10 0 

... 0 0-1 {(fc.-iyO + fc.Me-* 

... o 0 {(P,, - 1V0»* ♦ 1» -1 


Conclusions 

The effect of periodic accelerations on the stability of s 
multilayered, two interface, unbounded fluid system basbeen 
studied using Floquet analysis. In addition to the wave quid* 
her, five parameters appear in the problem: Bo it Bo,, Fr, p^, 
and p,,. Fr is inversely proportional to the square of the forc- 
ing frequency, and Bojifio J is inversely proportional to 

Several trends were discerned in the parameter study. For 
fixed density ratios fa and p,„ as well as fixed fr, the range 
of unstable wave numbers increases as Bo^BoJ increase*. If 
it is only the parameter Fr that is varied, it is found that the 
range of unstable wave oumbers increases as Fr is decreased. 
(Note that the variation in Fr values is very Grafted.) Physical 
interpretation of these trends has been presented in the pre- 
ceding section. 

Although the comparison is not presented graphically, the 
multilayered fluid system was found to be. In general, more 
unstable than the one interface fluid configuration. That is, 
the range of unstable wave numbers is smaller in the one 
interface case. In particular, the greatest contrast was in the 


The form of matrix A is given by 

...0 0 0 0 in 0 -tf ir* 

...0 0 0 0 0 in -k k 

...pi 0 0 0 p2 0 h (53 

...0 p5 0 0 0 p« p7 in 

with 

oJ - (1 + Pn + Pm). «2 “ (1 - PnW ♦ Pa) 

«3 - (1 - *0*1 ♦ P»). «< - 0 ♦ Pn) 

«5 * (I ♦ Pw), pi - (W2)e-»«2 
02 » (FrfiBoje*-^ l/a4, 03 - (Mr* *2 

04 - {Frfl)e-'a2, 0$ - (Fr/2)oJ OWQINAL PAGE 0 

OF POO* QUALir 
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£6 = - (fr 30o J )L'al'a5, p7 = (in) a3 
PS = (Fr 2)a3 
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Abstract 


with the increasing opportunities for research in a 
microgravity environment, there arises a need for 
understanding fluid mechanics under such conditions. In 
particular, a number of material processing configurations 
involve fluid-fluid interfaces which may experience 
Rayleigh-Taylor type instabilities due to external forcing. 
Ultimately, the stability of a multi-layer fluid configuration 
in a microgravity environment subjected to periodic residual 
forcing, both periodic and non-periodic, is of interest. This 
forcing may result from various sources, including equipment 
vibrations and space station maneuvers. 

An initial parametric study was performed to investigate 
the effects of steady forcing due to gravity; in particular, 
the limit case as gravity goes to zero. Solution of the 
linearized Euler and continuity equations and imposition of 
appropriate interface conditions yields a dispersion relation. 
This was solved using a root finder routine from the IMSL 
library. The numerical results are presented graphically 
using SAS/Graph via CMS. 

The utilization of a veil-integrated system of 
computational resources is essential for this research. 


Introduction 


The presence of gravity on earth is so omnipresent a 
phenomena that its effects often are not realized. In 
particular, gravity-induced problems arise in manufacturing 
processes (ie. buoyancy-driven convection in liquids, 
contamination from vessels that contain samples, and induced 
stresses that cause defects in crystals). An idealized 
laboratory would be free of gravity "contamination". 


(In Proceedings of the 1990 WVNET User Conference 
WV Network for Educational Telecomputing 
Morgantown, WV ) 
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The effect of gravity has been greatly reduced in 
low-gravity aircraft flights and drop tubes which provide 
short periods of microgravity, sufficient for some research, 
but certainly too brief for most processing experiments. The 
advent of extended spaceflight has dramatically increased the 
opportunities for long-duration research and development in a 
microgravity environment. Space manufacturing can eliminate 
the gravity-induced problems which are experienced in a 
terrestrial environment. 

The growth of crystals as electronic materials has not 
reached theoretical performance limits due to defects caused 
in part by gravity. During the spacelab missions, scientists 
were able to monitor crystal growth through each stage of its 
formation. In earth-grown crystals, it can be observed where 
the seed crystal stops and the new growth begins. The 
introduction of such a defect was not detected in space due to 
the lack of gravity- induced convection. 

The absence of convection is also pertinent to 
metallurgical manufacturing. A microgravity environment 
provides greater understanding of hov liquified metals diffuse 
through each other prior to solidification. Such knowledge is 
important for the production of improved and novel alloys. 

Containerless processing makes possible the production of 
much improved glasses and ceramics. In such a process, the 
sample is suspended and manipulated by acoustical and 
electromagnetic forces without the contamination of a 
container* Large samples can only be dealt with in a 
microgravity environment. 

Biological processing also benefits from space. Large, 
pure crystals allow analysis of many unknown protein 
structures which are essential to the design of new and 
Improved drugs. There is also effort towards the separation 
and purification of biological substances for pharmaceutical 
purposes* 9 * 

In the absence of gravity, fluid behavior which night 
normally be hidden by gravity-driven flows in a terrestrial 
environment, can be observed and analyzed. Drop and fluid 
column dynamics in microgravity permit experimentation of 
basic fluid physics theories* 

The environment of board space shuttle is not strictly a 
microgravity environment. Rather, residual accelerations 
exist which could affect any ongoing materials science or 
space processing experiments. A recent summary 12 indicates 
that accelerations Include those in the frequency range of one 
to ten herts, with acceleration levels ranging from 
10 *gwu to 10°*g«*rti>. In addition to periodic forcing 
(g-jitter)» residual accelerations may be of impulse type, due 
to such causes as station-keeping maneuevers and astronaut 
motion . 1 
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Most processes involve fluid dynamics, and in particular, 
fluid interfaces. This study does not investigate a specific 
process per se, but instead considers the stability of fluid 
interfaces. 

Previous work on fluid interfaces in a nicrogravity has 
focussed predominantly of the application of fuel slosh in 
tanks. Most recently, this has included work due to Hung et 
a l, which considered g-jitter in a slosh tank. A brief 
review of earlier work, as well as an extension of the 
efforts, was given by Gu. These investigations all involved 
liquid in a container of specified shape with a free surface. 

The stability of a single planar free surface subject to 
periodic forcing in the direction perpendicular to the 
interface has been investigated 3 ’ . Both studies were done in 
a l-g ambient environment and required the use of a container. 
In the work of Benjamin and Ursell , the container was 
cylindrical in shape. The analysis led to a Mathieu equation 
which governed the time-dependent amplitude of the 
disturbance. They were able to make statements concerning the 
interface stability based upon known mathematical properties 
of Mathieu equations. The case of a rectangular container has 
been addressed recently bu Gu et al , and the results extended 
into the nonlinear regime. Both of these investigations 
utilized an inviscid analysis. 

Viscous effects of the stability properties has been 
investigated recently in idealized infinite or semi-infinite 
configurations which had one fluid-fluid inter f aces. *'* The 
forcing was periodic and directed perpendicular to the 
interface. The work of Jacqmin and 0uval 6 assumed a zero mean 
g-level and pertains to a microgravity environment. A Floquet 
analysis was applied to the fluid system for the case of 
sinusoidal forcing. 

A multi-layer fluid configuration analysis has been 
performed by Lyell and Roh. The investigation considers a 
two-interface configuration in which the middle layer of 
finite height is situated between two semi-infinite layers of 
fluid. The analysis is inviscid and incompressible. A zero 
mean g-level serves as the base state for the study. 
Two types of time-dependent forcing are investigated, each 
simulating real microgravity environment accelerations 
(namely, g-jitter (periodic) and short-duration impulse 
(non-periodic) ) . 

Prior to investigating the time-dependent forcing on the 
multi-layer configuration, an analysis was performed for the 
case of a constant gravitational field. This particular 
investigation is the subject of this paper. 
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Problem Description and Equation Development 


Fluid configuration stability in the presence of constant 
acceleration fields is investigated. The limit cases of 
l*g«*rth and 0*g**rtn are studied, as well as various 
intermediate values. Lunar gravity can be approximated by 
0.16*g«»rth. The 0*g..rth state will ultimately serve as a 
basis for investigating the effects of residual accelerations. 

The configuration to be considered is comprised of three 
horizontal fluid layers. No rigid boundaries are present. 
The layers extend to infinity in the horizontal direction. 
The top and bottom layers are considered to be semi-infinite 
in nature, while the middle layer has a finite height. The 
geometry of the figure is given by Figure 1. 

The three fluids are considered quiescent; their mean 
velocities respectively equal zero. The fluids are immiscible 
and will be taken as inviscid and incompressible. Surface 
tension is a property of the interfaces. A normal mode 
approach to the perturbation is assumed. In a normal mode 
approach, the disturbance (or perturbation) is assumed to be 
wavelike. If the wave grows, the fluid system is said to be 
unstable to the disturbance (perturbation) . If the wave 

decays, the fluid system is stable to the perturbation. This 
is mathematically represented as follows: 

, - • ( * ike « t) • (ke i t) (l) 

where n ■ interface shape 
€ ■ amplitude 
k - wave number 

c K > real component of propagation speed 
c f * imaginary component of prop, speed 

A positive value of c ( will cause growth of the perturbation; 
and hence, the fluid system will be unstable. 


Two cases are investigated: 

CASE 1: air/silicone oil/water (region 2/region 1/region 3) 
(stable configuration in terrestrial environment) 

CASE 2: air/water/air ■ 

(unstable configuration in terrestrial environment) 

The parameters to be varied include height of middle slab 
(h) , wave number of the perturbation (k) , and the value of the 
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constant gravitational acceleration (g) . The propagation 
speed of the perturbations can be calculated for different 
parameter conditions. A positive value of the iaaglnary 
component of the propagation speed will indicate an 
instability on the fluid system. 

The governing equations of incompressible fluid mechanics 
are the continuity equation and the conservation of momentum 
equation. The analysis is inviscid and linear. Linearization 
is done about a state of zero mean motion. The following 
linearized equations govern fluid behavior. 

7*u = 0 (2) 


<3u 

p — = -Vp + pge 

at 


(3) 


where u » velocity field 
p =» pressure 


It is straightforward to solve these equations for both 
the velocity field and pressure. (A velocity potential 
formulation is utilized in the solution.) 

The dispersion relation is obtained by applying three 
boundary conditions at each interface. These three conditions 
are: (1) the kinematic boundary condition, which states that 
a particle of fluid which is at some point on the surface will 
always remain on that surface, (2) the matching of the normal 
component of the velocity across each interface, and (3) the 
normal force balance across each interface. 

The boundary conditions are applied to the governing 
equations, and via manipulation, a dispersion relation is 
obtained. It is given as follows: 


(P 1 +P 3 ) (P x +P 2 )e 


2kh 


(P 2 -Pl) (P^Pj) 


p 

{<P t +P,)( f(Pj-P l )-T 2 )c) + ( f(P 1 -P 3 )-r,k)} 


,21th 


+ (P 2 -Pj)( flPj-Pj)* r 3 *) + (Pi'PjK l(P 2 -P x )-» a k) 


C 


2 


019 



« 0 


( 4 ) 


► 


2kh 


( ?^2“ p l ) ” T 2 k H ^ (P 3 -Pi> -y 3 )c) 


Thus, this propagation speed, c, is the solution of a 
fourth order polynomial. It is also the eigenvalue of the 
stability problem which is posed by equations (2) and (3) 
together with the associated interface conditions. The 
propagation speed is a complex number, c - c ♦ ic . 

An imaginary component equal to zero implies a neutral 
disturbance. If Its value is less than zero, the exponential 
term decays in time, and the system remains stable. However, 
if this imaginary component, c j( is positive, the exponential 

term grows in time, resulting in an instability. 

A limit case which is analytically tractable can be 
obtained from the full dispersion relation for the special 
case in which the ratios of the top and bottom densities to 
the middle density are negligibly small. For such cases, the 
configuration will remain stable if the following inequality 
holds true: 



* 



( 5 ) 


The scope of this study is to analyze the four 
previously stated cases under various parameter conditions. 
That is, by allowing the parameters to vary over a specified 
range, the roots of the dispersion relation can be calculated 
and, hence, Interface stability can be determined. The 
parameters that are considered are the height of middle 
layer, the wave number, and the value of the gravitational 
constant. For our ultimate purposes, we are most interested in 
the case in which the time-independent gravitational body 
force is zero. 

For a more detailed description of the equation 
development, see Roh. 
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Utilization of WVNET Resources 


The numerical solution and graphical representation of 
this analytical problem requires a well-integrated host of 
computer resources. The CMS system was accessed via WHET on 
a VT320 terminal. A remote sight at the Engineering Sciences 
building was used. 

The dispersion relation was solved numerically using the 
DZPORC routine of the IMSL library. The DZPOR£ routine makes 
use of the Jenkins-Traub three-stage algorithm 7 , in which the 
roots are computed one at a time for real roots and two at a 
time for complex conjugate pairs. As the roots are found, 
the real root or quadratic factor is removed by polynomial 
deflation. 

Two options were explored for graphing the results. 
Initially the data was downloaded to a diskette via Kermit, 
which in turn was plotted using Lotus/123 graphics package on 
a Zenith DS computer. While the output was satisfactory, it 
was inconvenient and time-consuming to change terminal sites. 

The second, and preferred, option was to access CMS 
directly through a WVMET line connected to a Macintosh II PC. 
This was accomplished via VersaTera and VersaTerm Pro 
communications. The program calling IMSL routines was run as 
simply as with a VT320. The data was then transferred to a 
SAS/Graph routine emulating TEK4014 device, which presented 
the results graphically. The plots were converted to MacDrav 
files from which hardcopies were obtained. The advantage to 
this option is the one-terminal site capabilities. 

The numerical results for the time-dependent forcing also 
accessed several routines from the IMSL library. One 
solution, in the case of periodic forcing, involved the 
eigenvalues of a very large complex matrice system. An 
enormous amount of storage space was required for computation. 
Temporary disks had to be accessed to provide the necessary 
space. Details can be found in Roh. 1 


Results and Conclusions 


The fourth order polynomial (in c) has four roots. 
Because of the nature of the dispersion relation, the roots 
were generated in pairs. That is, for any given solution, 
there exist two pairs of roots, where each pair consists of 
the positive and negative values of a number. Physically, for 
real roots, this means the perturbation may propagate in 
either the positive or negative direction. For imaginary 
roots, it implies an instability will occur since these roots 


D21 




occur in complex conjugate pairs. If all the roots are real, 
the system will be stable. 

Figure 2 shows the four roots of the dispersion relation 
for Case 1 (air/silicone oil/water). The roots are plotted 
over a range of gravity ratios (g/g*«rth) from 0 to 1.0. As 
might be expected, this case is stable for all parametric 
conditions. (Note that the non-zero roots are exclusively 
real.) A less dense fluid above a more dense fluid is stable 
to small perturbations in the presence of a constant 
gravitational level. 

Figure 3 show the dispersion solution for Case 2. For 
certain gravity ratios, the presence of a positive imaginary 
root is revealed, which implies an unstable configuration. 
This behavior is expected since the configuration has a more 
dense fluid above a less dense fluid. 

Since an instability depends solely on the presence of 
positive imaginary roots, the subsequent figures will display 
these particular roots exclusively. 

The effect of wave number, (k), on configuration stability 
is elucidated in Figure 4. As k values increase, the 
configuration becomes more stable. Thus, since k is inversely 
proportional to wavelength, the configuration is unstable to 
long wavelength perturbations. The restoring force required 
to maintain stability is greater in the long wavelength 
regime. Note that all cases are stable at 0*g««nh. The 
results of Case 1 do not appear since the configuration is 
stable for all parameter space. 

From Figure 5, it is tempting to conclude that the middle 
slab height, (h) , has no effect on the stability of the 
configuration. This conclusion is valid for values of h which 
are considerably larger than the wavelength (inversely 
proportional to k) . However, as h approaches the order of the 
wavelength, the effect of the slab height becomes apparent 
(see Figure 6) . The fastest growing instabilities are 
associated with configurations with the smallest middle layer 
heights. 

The limit case (eq. 5) simulates a liquid layer situated 
between two layers of a gas, and its accuracy can be verified 
by comparing it to Figure 4. According to (eq. 5), for Case 2 
(air/water/air), and h»lcm, the instability should originate 
at g/g*«rth»0.073 for k*l, at g/g«*rth« 0.294 for k«2, and at 
g/g««rth*0.661 for k»3. The results from Figure 4 confirm 
these values. 

It is seen that in the case of zero gravity, each 
configuration remains stable. Although we might expect 
Rayleigh-Taylor type instabilities for Case 2, there is no 
body force which would drive the density difference;hence, the 
system will remain stable. 

This zero gravity state has been taken as the base state 
for the time-dependent studies. As previously stated, the 
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presence of g-jitter or manuever type short-duration impulses 
will play an important role In the stability of 
multi-interface configurations (see Roh) . 
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Copy of Course Syllabus: MAE 399 



MAE 399: FREE SURFACE AND INTERFACE FLUID MECHANICS 


L ^A SQN EQR i^ IffiSLfQRJ^EjgS 


The student (Ms. K. Perkins) re quires knowledge of free surface fluid mechanics. 
Moreover, it is necessary that the student learn this material during Fall Semester 1993. 

This course has information that the student will need to know even before being told of her 
thesis problem. No other MAE class covers the material which will be taught in this 
directed study class. As the material has a large mathematical component as well as 
demanding interdisciplinary fluid physics, it is felt that a structured class is necessary. 


The student will be working on a thesis topic in the area of free surface flows in a 
microgravity environment. The topics covered in this class will provide some (not all!) of 
the necessary background knowledge. 


2. NARRATIVE DESCR 
PREREQUISITES 


limW 


Co-current registration in MAE 411, 
Background of MAE 104 and 
Co-current registration in MATH 317. 


DESCRIPTION OF COURSE 

This is a course in free surface fluid mechanics. It includes both invisicid and viscous 
flows. The treatment will be equation-based, with both analytical formulations as well as 
pertinent numerical method approaches to be covered. The conservation of mass, momentum 
and energy equations will underline the work, of course. In the case of the inviscid 
formulations, the topics covered will involve traveling surface waves. If the geometry is 
finite (and the fluid is containerized), the focus will be on standing waves. For the viscous 
fluid formulation, many practical problems will involve a thermal field and/or a 
concentration gradient. The analytical formulation becomes much more complex. 

Marangoni flows become very important, and will be investigated. 

Some of the physico-chemical aspects of the free surfaces will be introduced, beyond 
the basic concept of surface tension. Concepts such as excess surfactants, absorption and 
des^tion, and reaction at the interface will be introduced. This will necessitate the 
introduction of the diffusion equation (for mass) as another pertinent governing equation. 

The description of the free surface as an interface will be introduced, along with such 
concepts as surface viscosity, etc. 

The application of free surface flows in microgravity will be covered, also. This will 
include liquid column (float zone) applications, planar configurations involving Marangoni 
flows, and a brief introduction to drop dynamics. 

Finally, the current state of the art in the computational approaches to resolving wave 
motion on the interface/free surface and in the determination of bulks flows with a free 
surface will be elucidated. 
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3. TOPICAL COVERAGE 


l. Free Surface Inviscid Flows 

A. Governing Equations -Inviscid Flows 
Boundary/Interface Conditions: General 

B. Wave Motion - Linearized Problem 
Formulations/Linearizadon Scheme 
Dispersion Relation - Travelling Waves 
Standing Waves in Finite Containers 
Waves in Forced Containers (Intro) 

C. Nonlinear Wave Motion, Solutions 

n. Free Surface Flows/Viscous Effects 

A. Extension of Governing Equations 

Discussion of Restriction to a Newtonian Fluid and Alternative 

B. Boundary/Interface Conditions 

C. Discussion of Physical Examples 

D. Introduction to Computational Methods 

m. Free Surface Physico-Chemical Hydrodynamics 

A. Properties Beyond "Surface Tension" 

B. Introduces posibilty of longitudinal waves as well as transverse waves 

IV. Free Surface Flows/Thermal Fields 

A. Effect of Thermal Gradient on Surface Tension 

B. Marangoni Flows 

C. Physical Examples, Introduction to Microgravity Applications 

V. Free Surface Flows/Concentration Gradients 

A. Surfactants Confined to Interface 

B. Surfactants Absorbing/Desolving in the Bulk 

C. Effect of Concentration Gradient on Surface Tension 

D. Application/Microgravity Examples 

VI Additional Computational Methods 

A. Emphasis on Algorithms, Approximations Made (eg. How much does the 
interface deform?) 
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Topic I: 

a) I.G. Cupkie 

Fundamental Dynamics of Fluids 

McGraw-Hill 

Chapter 6 

b) G.B. Whitman 

Linear and Nonlinear Waves 
Wiley 

Chapter 11, Selected Topics of Chapter 12, Chapter 13, Sections 1-9,10-12 

Topic II: 

a) Dr. Lyell’s Notes 

b) Computational Fluid Mechanics 
Chow 

McGraw-Hill 

Sections 2.9, 2.10, 2.11, 3.6, 3.7, 3.8 

c) Papers from the Literature including: 

H. Bauer’s on Natrual Oscillations of Viscoelastic Fluid Drops and/or Cylinders 

Topic HI: 

a) B. Levich 

Physico-Chemical Hydrodynamics 

Article in the 1963 Annual Reviews of Fluid Mechanics 

b) T. Sorensen (Editor) 

Dynamics of Fluid Interfaces 

Chapters (actually tutorial papers) by Hennenberg and Sanfeld 

Topic IV: 

a) Dr. Lyell’s Notes t 

b) Stemling and Screen’s Gassical Paper 

c) Recent Papers from the Literature which Involves Microgravity Formulations (and 
thus neglect the complicating bouyancy force) including: 

* Doi & Koster 
Phys. Fluids, 1993 

* Myra Carpenter! 

MS Problem Report/WVU 

d) Scientific Foundations of Space Manufacturing 
Avdeyevsky 
MIR Publishers 

Topic V: 

a) Scientific Foundations of Space Manufacturing 
Avdeyevsky 
MIR Publishers 

b) JFM Papers by 
-Sen & Davis 
-Meiburg & Homsy 

on Marangoni Flows in Slots with Contamination 

TejatSL* fro™, H*. UfanJun l PJujvw* 

Ack'l! Q. Ot iv. tluJL, 0. 

WA. firMJC d E3 



4, GRADING 


HOMEWORK 

40% 

TEST 

25% 

COMPUTER PROJECT 

10% 

REPORT 

25% 



